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FOREWORD 


PART I 


7b la Is a combined report for the two projects at the Massachusetts Institute of 
: - or r.ology Wf '.ch ere sponsored uy the Office of Naval Research under Contract N5orj60. 


: rojec* on Machine Methods of Computation and Numerical Analysis 

Th‘.s roject la ar. outgrowth of the activities of the Institute Committee on Machine 
K. tnida :f Conpu a : .'.or,, ea^atilshed In November L950. The purpose of the Project la (l) to 
• y-atT the efforts or a., the departments and groups at M.I.T. who are wording with modern 
nj 111 rtg machines and their applteat tors, and (2) to train men In the uae of these machines 

for computation and numerical analysis. 

People from several departments of the Institute are taking part In the project. 
r» the Appendix will bo found a list of the personnel active In this program. 


Project Whirlwind 

This Project makes use of the facilities of the Digital Computer Laboratory. The 
prif'j’pal objective of the Project is the application of an electronic digital computer of 
.arge capacity <ir.d very high speed (Whirlwind I) to problems In mathematics, science, engi¬ 
neering, simulation, and control. 


The Whirlwind I Computer 

Whirlwind I Is of the high-speed electronic digital type, In which quantities are 
.-.•presented as discrete numbers, and complex problems are solved by the repeated uae of 
fundamental arithmetic and logical (t.e., control or selection) operations. Computations 
are executed u;/ fractional-microsecond pulses In electronic circuits, of which the princi¬ 
pal ones are (i) the flip-flop, a circuit containing two vacuum tubes so connected that one 
tube or the other Is conducting, but not both; ( 2 ) the gate or coincidence circuit; ( 3 ) the 
magnetic-core memory, in which binary digits are stored as one of two directions of magnetic 
flux within ferro-magnetlc cores. 

Whirlwind I uaea numbers of lb binary digits (equivalent to about 5 decimal digits). 
This length was selected to limit the machine to a practical size, but It permits the computa- 
' ion of mar.., simulation problems. Calculations requiring greater number length are handled by 
•.fie ise of multiple-length numbers. Rapid-access magnetic-core memory has a capacity of 
32,768 binary digits. Present speed of the computer is ^0,000 single-address operations per 
second, equivalent to about 20,000 multiplications per second. 


Machine Methods of Computation and Numerical Analysis 
1. GENERAL COMMENTS 


Two new studies originating from the mathematics members are reported here. Both 
concern mathematical techniques in the study of deflections or shallow shells. In one case 
the convergence of iteration procedures to the solution of a differential equation is the 
subject. In the other case the objective la to obtain integral representations, asymptotic 
behavior and recursion relations for certain functions In order to solve boundary value 
problems• 

A project that has been completed is the evaluation of multi-center integrals. This 
work that will be reported fully elsewhere la used In a calculation of energy bunds for a 
graphite crystal model. 

In the way of theory of programming and general use of a digital computer, a complete 
report on the programming algebra introduced last time Is given. Though this algebra la a by¬ 
product of a general approach to minimizing computing machine time, it serves some practical 
purposes Independently of the general theory. The place of this algebra In reference to the 
mathematical definitions of computing machines la also discussed. 

The Monte Carlo technique for Implying diffusion behavior from the study of individual 
neutrons has been tested throughout a number of component subroutines and Is close to being used 
for reactor studies. A detailed outline of the logical design for this Monte Carlo method la 
Included In these reports. 

Six other projects that have been reported before In their early stages are progress¬ 
ing. These will bring results In the fields of physics, geology, civil and mechanical engineer¬ 
ing and concern Coulomb wave functions, earth resistivity, dynamic loading on buildings, crack 
growth In bars, and vapor condensation. 
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2. ORADUATE SCHOOL RESEARCH 


2.1 Index to Report# 


Title Pe«e 

The Beale Problem of Numerical Analysis Expressed In the Language 

of Computing Machines 7 

Finite Bending of Thin, Shallow Spherical Shells 11 

Properties of the Solution of a Differential Equation Arising In 
the Theory of Linear Elastic Deformation of Shallow Shells or 

Revolution 3 3 

Asymptotic Solution of a Differential Equation lil 

Evaluation of Multi-center Integrals 15 

Coulomb Wave Functions 15 

First Approximation Solution on Ore Body 16 

An Application of Monte Carlo Methods to Neutron Diffusion 20 

Response of a Single Story Reinforced Concrete Building to 

Dynamic Blast Loading 22 

Response of a Five-Story Frame Building to Dynamic Loading 22 

Growth of Patigue Cracks 22 

Condensation of Vapor in a Vertical Tuba 23 


2.2 Progreas Reports 

THE BASIC PROBLEM OF NUKKHICAL ANALYSIS EXPRESSED IN THE LAN0UAOE OK Mt'UTlNO Mm'HINK 


The first two reports on this subject [l,2] led to a prig ramming aigebr and in trie 
second of these reports we stated that the algebra was complete except for dec is. n making. 

The remaining essential point we Intend to discuss this time. 

DECISION MAKING 

It la Largely the decision making that gives to computer programming a unique aspect 
not known In other branches of mathematics. In 19*7 von Neumann and Ooldstlne I'l assoc idled 
the words ''dynamic'' and "macroscopic" to this part of the programming. In their words, “Since 
coding la not a static process of translation, but rather the technique of providing a dynamic 
background to control the automatic evolution of a meaning, It has to be viewed as a l glca 1 
prublem and one that represents a new branch of Torma1 logics." Using this term 1 no logy we 
can aay that the past progress report was concerned largely with the static translation «*i a 
problem. The presence of the modification instruction T mnD and the po 3 lDracket with a rained 
Integer ( ) n served only to condense a program from a direct transliteration or a problem t 
an array of subroutines. The machine made no choices; all decisions about the values ul the 
subscripts on T and the raised Integer after ( ) had to be made at the time of the progromm'ng. 
The conventions of the algebra are constructed, however, to permit these subscripts ami r .! **J 
integers to be equated to the outputs ol' other programs: that Is. Instructions in one program 
can be modified by the outputs of other programs and this Is the first step toward containing 
a decision logic within the algebra of programming itself. 

In the present report we will generalize the possibility of modifying instructions. 

At the same time we will Introduce the basic decision Instruction. From the mathematical 

point of view the decision Instruction plays the part of a function that Is 1 over a set and 

zero elsewhere; it combines with ( ) n to play the part of a conditional transfer In conven¬ 
tional programming. 

First, let us review* the Instructions already available and at the same time use 
a more convenient notation. Throughout this report 0( ) will stand for the contents of the 

last register referred to in the sequence of Instructions contained in ( ), unless 0 is a 

vector in which case we will write 0 M . 

a: The identity. Replace contents of register a by Itself. 
a^ c : Add 0(b) ♦ 0(c). Store In a. 

a b-c : Subtract 0(b) - 0(c). Store in a. 
a b*c : 0( b ) • 0(c). Store In a. 

a b ^ c : Divide 0(b) / 0(c). Store In a. 

T mnp : If P receded b y a b+c* a b-c* etc ** re P lfl ce b by *►+«". c by c*n, a by a*p. 

( )nj Repeat bracketed sequence n times. 

The generalization or the modification Instruction operating on la a set 

of equations as follows: 

*T+z 

x - o( j ) 

Y « 0( „ ) : Perform the three sequences of Instructions contained 

z 0( Z ) ln brackcts 1*2,3 ln tbat order and than vrtorm Xy +Jf . 

A similar definition holds Tor the modification of other Instructions. (Techn1ca1ly, this Is 
not a strict generalization of the simpler mod 1 Pleat ion instruction, because '.n the simpler 
case the Instruction a. was followed by modification rather than preceded by It. The 
difference is not lmportSnt.) 

The generalization of the repeat Instruction is a set of equations as follows: 


* Note: Gome notat Iona 1 errora were present In the definitions given last time. 
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0 ( 0 ) 


( , >" 

H • 0( 2 ) 


: Perform bracketed sequence Z. Then carry out 
Instruction ( j ) . 


The decision Instruction is defined as follows. We will assume throughout that 

0 , 0 ( 1 ) - 1 . 


(*0+0 lf 0(b) ‘ 0(c) 
b * c [a^ IT o(b) »0(e) 

We will also uee the sywbol 0 aa the Inverse of 0: 


mappings (vector and single valued, respectively) denned on this space. Por each point 
(p,i) In PXI, 0„(p,l) and T(p,l) represent the output and time. Thin definition differs from 
that given by Turing In a number or respects. First of all, each point (p,l) In CXI 

orresponda to a Turing machine for which 0 H (p,l) - {<*,.11,,.. .<i„) Is a vector representation 
the associated Turing computable function. (Strictly speaking, since for each (p,l) there 
la a finite collection or numbers Q, ,d 2 ... .ttw defined only Tor the finite sequence of Integers 
1,2,... M, each computable function la partia l. The fact that H, are rational whereas the out¬ 
puts of a Turing machine are Integers Is Inessential.) Second, since Tor each point (p,l). 
n„(p,i) and T(p,l) are denned (and thus finite) the machine R can never bee roe Involved In 
an Infinite loop. In logical terminology, the machine M can never become Involved In an un- 
deetdable propoaltlon. Thus In defining M we exclude from the atart tne possibility of 
handling problems that esnnot be resolved In a finite amount of time. 


0(a) • register containing a. 


Thus, lf 1 • 0(a), > • 0(X). 

A final Instruction Introduced to complete the algebra Is 

a : Take the coefficient of 2 " 0 ( c ' in the binary representation 

b.2“ c of Q( b) and store It In a. 


The generalized modification Instructions (Instructions made up ..f variable registers) 
and the generalized repeat Instructions (brackets with variable raised Integers) may appear In 
rather complicated combinations. For Instance, there may be several bracket:-, within one another 
all accompanied by different variable raised Integers and within all nr some of these brackets 
there may appear Instruct Iona with variable registers or variable contents of registers. Un¬ 
less a convention Is established. It will be ambiguous with respect to w-uch brackets the In¬ 
structions are variable and with respect to which they are constant. For example the program 


t (X 3v .v> N l 


n 


may mean that the X,,X_,X., are to be evaluated each of the M timeii the outer or. c<*;t U. entered 
aa well aa each of ihe^N times the Inner bracket Is entered. Alternatively, It may mean that 
the compound Instruction 


(X. 

V X 2 


) 


N 


appearing Inside the Instruction ( ) 


n 


Is to be considered constant and 


There Is another Important distinction. The basic problem that we have studies In 
these reports Is that of finding the most efficient program for coming within a certain dis¬ 
tance ol a dealred output. In formulating thla problem we have hud to assume that at least 
one program la known (efficient or Inefficient) whose output characteristics are desirable. 

In other words, we have bypassed the problem of greatest Interest to the logicians; that is, 
the problem of what are the computable functions, or what are the problems that can be pro 

grammed and solved on a machine In a Tlnlte time. Before we can approach a solution to either 

of these problems In the context of the present definition of M, it lo necessary t have a 
realization of N. Such a realization la, In fact, supplied by the algebra of instructions 
Just presented. It Is clear that the conventions of the algebra yield all finite sequences, 

1, of Instructions and thla collection can be called I. The convention that all elements or 
the sequence p are to be the successive contents of the register numbers 3,k,5... appearing 

in the algebra completes the realization of PXI. If the components of the vector Ogfp.l) 

a re defined as the rinal contents of any finite col lection of registers depending on (p,l), then 
0g(p,l) 1* determined by (p,l). T Is a free variable and may be specified Independently. The 

algebra together with a specification of T thus brings ub a realization of M. At the same 

time we have a tool for programming that reflects the mathematical structure of the problems 
being solved. 

A complete solution to the problem of minimizing performance time awaits a fuller 
investigation of the laws of the algebra, of which one or two were discussed last time. The 
first solution to this problem will come when the rules of the algebra can be used to explore 
completely the region In PXI that la algebraically equivalent to an arbitrary point (p 0 »1 0 )» 
From then on. It will be necessary to select out of each algebraically equivalent region 

the point or points that minimize T. A solution to the logically Interesting problem of 

characterizing the class of problems computable on H also await further Investigation of the 
Instructions that appear In the algebra. 


EXAMPLES 


Is variable. Thus, the convention will be as follows: Beneath the equation for computing 
any raised Integer (and lf the integer la constant, we enter the equation as N - N) will be 
listed in order all of the equations for registers and contents of registers that are variable 
with respect to the corresponding bracket and raised Integer. (See examples.) If a variable 
with respect to one bracket is constant with respect to a second, the equation for Its evalua¬ 
tion will not appear In as oclatlon with the raised Integer of the second bracket. The order 
of the equations that do appear In conjunction with a bracket should be consistent with the 
order of appearance of the variables within the corresponding bracket. 

A program written In this algebra will appear as a sequence of Instructions together 
with aa many sets of equations for evaluating the variables as there are brackets. By follow¬ 
ing the above rules the temporal Sequence for performing Instructions will not be lost, even 
though, for convenience, certain parts of the program are written as separate equations. The 
rules for combining the primary sequence of Instructions with the accompanying equations and 
the converting to the instructions of an ordinary comoutlng machine will be discussed after 
some examples. 


DISCUSSION 

Some remarks on the generality or this algebra and Its relation to the definitions 
ol computing machines given in these reports [1,2 ] und elsewhere will be appropriate. The 
mathematical definition ^f a computing machine M that we have used Is M ~(PXI,O n ,T) where 
PXI Is a space of finite sequences of parameters and instructions, and 0^ and T are real 

4 Note: We have not previously specified that the sequences In PXI are Tlnlte, but this 
restriction will be assumed from now on. 
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As examples of how the algebra might be useful as a programming tool. Irrespective 
of tne general problems discussed above, we consider the following that arose in the electron- 
photon cascade problem [5] 

F.xample 1. It ^s required to gene rate the array of numb ers II from the ar ray I and pu t them 
In the registers d laplayc? Tn t he array~TII . The program shoul d be suTtahl e for sl ml lar 


— 

- *- - 



m 
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let 
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0 

O 

o 

• 22 

nt 

in 

167 

It 

Is seen that the array 

II 

Is generated 

by addition In an obvious «ay 


lem Involves one basic addition Instruction 


'V X 2 

that will be used a variable N number of times for each column. The cycle corresponding to 
each column must be repeated a fixed number M times, where H+l Is the number of columns. There 
will be a special program for changing the alia of N. 

There will be two aeta of equations for computing the variable values of X. , X-, X. because 
these variables ohange between rows and between columns according to different rbles. As It 
hardens, the variation between rows affects all three and the variation between columns affects 
only one of the X^ f X 2 , X^ . Thus there are four auxiliary sub- 
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routln.. to th. program or brook.t for ro... on. to oo^uf H ud thr.. to co-put. th. T.rl.bl. 
r.Rl.t«r« X,, * 5 . x.. Thore Is one auxiliary subroutine to the 

bracket Tor colulna. We .111 use the ay-bole al.a2.a3.** ae floating aOdreeses (that can 
later be given determined values) Tor the registers containing th.currentvalueaofthe 
Tour variables. The Initial values of the four variables must be determined by the boundary 
conditions of the problem. 


Tne program can be written aa follows: 

Y 


f (X, )’ 

1 5 X l+ X„ 




H 

0 ( »*„ 


°<* l al + l) 

0(a \.4-l> 

0(a2 a2-l> 


H - ° (a ’.4el> 

In reference to arrays I. II, III, the initial contents or registers 117, ai, a2 and ak 
should be i, 2, 121, 118. The initial contents or a 3 can be anything, M should be set M - 3. 

Example II. Let 16 binary fractions or lb digit accuracy each be stored In registers 101,... 
lib. let us examine tfie first H, dlgfEs ^o the~rlght ofthe decimal ip the fj.rat number. In 
FHe second number, etc., recording tne ~ ToEa l number of^Tlglis 1 round . 5t°g a fter examining 
aTT lb numbers or after examining enough numbera to make the total number of 1 *a greater than 
>r equal to 2. 

The basic program here Is to search and store. There will be a bracket for selecting 
a new digit, a bracket of selecting a new number, and a bracket for stopping the process after 
sixteen numbers have been examined. The digits selected will be variables and the registers 
containing the numbers will be variables. Since the digit selected must be re-set to zero 
(a re-set parameter) after a new number has been chosen, the digit number appears as a variable 
with respect to two brackets. One of the raised integers must, of course, depend on the rela¬ 
tive number of digits 1 found. The program Is as follows: 

IL 

; r (ai _ 0(x , .5 ) l ] p 

l I X r 2 0 'V a>al 1 1 


*2 * "2 


°<* 6 a4,aj> 

Q(a5 o+o) 




To aatLafy the boundary conditions, the Initial contents of registers a2,a3,a4 should be 100, 
0,2. The initial contents of al,a5,ab may be arbitrary Mg should be set Mg - 16. The 
result will be stored In a3. 


CONVERSION OF THE ALGEBRA 


ning 


The first Btep toward converting the algebra to any traditional form of program 
Is the elimination of all the variables except the raised Integers. This Is done by (1) 
insert lng. all the programs for computing the variables In order of appearance of the variable 
Into the left-hand side of the corresponding bracket and (2) replacing all the variables by 
the symbol 0(a), where a Is the register containing the variable. At this stage In Example 1 
the program will appear as follows: 

N.M 


[ a4 a4-l (al, a4-l aZ a2-l aJ a4+l 0(8,1 >0(a2)+0(a3) 1 l 


M - M, 


0(al aHl> 


The second atep Is the elimination of the variable ralaed integers. Thlo la done 
by (1) Inserting the program that computes the variable Immediately before the corresponding 
bracket and (2) replacing the variable by 0(a), where a Is the register contalnlng the variable. 
In the above example we would then Insert al ... Immediately arter a* , , and before the Inner 
bracket. N would be replaced by 0(al). al+i a *‘ l 

This much done, the third atep, the translation of the linear array of Instructions 
of the algebra Into a linear array of given Instructions remains. The details of this step 

10 


graduate school research 

will not be dlecusaed because each machine has Its own apeclal Instruction. However, a row 
remarks may be helpful. Instructions Involving the symbol 0(a) m.at be translated with the 
a id of a mod 1 flea t Ion Instruction such as “tranarer tne contents of register a to the address 
part or Instruction — —. Both the decision Instruction and the repeat Instruction will 

Involve conditional and unconditional transfers of control. The repeat Instruction win 
certainly Involve a counter of some sort. The Instruction Tor examining a digit will usually 
be translated with the aid of a cycling Instruction. 

18 clear that one to one correspondences can be made between the Instructions of 
this algebra and sequences or Instructions Tor any high-speed digital machine. Once these 
correspondences have been established, the conversion of the algebra to any machine language 
la purely mechanical. 

Bayard Rankin 
Raymond R. Stora 

References: 
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Proc . Lond . Hath . Soc .. Series 2 , 24 (1936) p. 230-265. 

[5] Raymond Stora, Same as [l] , Report Ho. 18, December 15 (1955), p. 12. 

PIHITE BENDING OF THIN, SHALLOW SPHERICAL SHELLS 

1. The problem under consideration Is that of amall finite deflections of thin 
snallow spherical shells of constant thickness h subjected to a uniform radial load p. The 
differential equations of this problem have been derived by F.. Relssner tl). They can be 
written In the following form [2] : 


(1) H[F,0] - F” + i P> - ij P + ^ 2 a . Zx - 

H (P, 0] 


0" + 


1 

? 


r *♦ rh 


where y Is proportional to the load and>* la a geometrical parameter which is zero for flat 
plates. I^t ▼ be a stress function [1J and Q be the angular deflection of the shell, then 
P, G, y % and are defined as follows: 

.jV 

4Eh a_1 ‘ 

where 


9>»! . o.-V,. . 2 

7* Eh^ ' 4Eh^ ' “ 


a * radius of the middle surface of the shell 
2<* • opening angle of the shell 
E m Y oung's mod ulus of elasticity 

m Z - V\2(l - l?) 

V • Poisson's ratio. 

We shall consider two cases of boundary conditions [2]: 

(2) A simply supported edge: F'(l) -^P(l) • 0• (1) -*>G(l) • 0 

(3) A clamped edge: F * (1) - O'(l) -J>0(l) - 0. 

Observing that x - 0 Is a singular point of the system (1), we require that F and 0 be finite 
In both cases (2) and (3): 

II 
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(») p(o) c °° 


0 ( 0 ) ^ 


ye .tte.pt to find .elution, or equation. (1). (Z) and (A), or (1). (3) and (A) respectively, 
in form <>f power aerlea expansions: 

p . ^ P *" 0-21 0 *". 

|j»# n a • • 

Substituting these series Into the differential equations and equating the coefficients or 
like powers of *, one finds the following recursion formulas 


(3) Z, - £(Z -^ 2 °i ♦/l' 1 0 1 ) 

P n • ±1 <-/Vz * n - 5 ' 7 ' 9 . 

°n * </Vz ' i P /n-J-l> " - 5 ’ 5 ' 7 . 

F zk ’ °zk ' 0 k - 0 ' 1 ' 2 . 

If F. and 0. are determined, the remaining coefflclenta are given by these recursion formulas. 

Tne stresses and displacements may also be expressed In terms of P R and 0^. 

2. To obtain initial estimates for the coefficients Pj and 0^, we approximate the solutions 
P and 0 by 

P # - Ax ♦ A»x 5 ; 0* - Cx ♦ C»x 5 . 

In order to determine A, A’, C, and C 1 we shall require that P and 0 satisfy the differen¬ 
tial equations (l) in the mean;, 

j I 

/H[p', 0 *]dx - 0 and /H[P*,o'Jdx - 0 , 
z> O 

and also that P* and 0* satisfy the boundary conditions (Z), or (3) respectively. Using these 
four conditions to determine A and C and taking 

r. A, » C, we obtain an approximate solution F^ # G^°^ by 

means of the recursion formulas (5). In order to obtain successive improvements of the 
coefficients P and Q , the following procedure.has been proposed by R. Simons [2j . Denote 
the mth approx^matlon n of (P_,G ) by (p \ m J) and let 

(6) Pn (-i, n . Jw + iy, o n <-) - a n <") * Mn ( » 

P and are required to satisfy the recursion formulas (5) except for terms of 

higher order than one In the corrections APk'?J. sn<) AG*'®' Hence , one obtains recurrence,formulas 
for and which are linear ln^P^* and d0 k ( " ,, (k*l..n). Therefore A 1 and 

A0 can be written as linear combinations of^F.^ and doJ m . 

(7) 6F n W - W mi 

rtj 9) - . « AF <■> ♦ d <*> AO 
n n l n l 

where a b ^ I ”', c ^ and d are,given bv recursion formulas similar to ( 5 ). involving , 

only the known coefficients F and 0 (”). The second requirement on (P (m+1 J G_I ■"rl J) 

la that n n ' n n 

Zp(" Hl )x"and fo<“ + 1 )x n 

h h 

satisfy the boundary conditions (2), or (3) respectively. Using (7) this can be written 

( 8 ) AF 1 (b) Z('' ♦V)» n ( " ) ♦ fiO <">/;(« *i,)b n (B) - - T [n 

4 * n 1 h n * n 

AP 1 (B) Z( n -* , )e n ( " ) ♦ AQ l < " ) ^' < n -^) 8 n (B) - - £> -t') 0 n (B) 

it remains to solve the linear system ( 8 ) to obtain the (o>+lHh approximation P q ("t- 1 ) 

aa given by (6). Finally the coefficients p (m+1) . Q (mil) (1>1) oomput ad from thi recur- 

slon formulas ( 5 ). ' t 
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■. Computations are now under way on the Whirlwind I Computer for various values of the 
aramrters V and In particular we wlah to determine the range 1 parameters y andM 2 

or which the above Iteration procedure will converge to the solution of ihe differential 
equations satisfying the appropriate boundary conditions. We have obtained solutions fur 
o.o I live and negative values of V In the range lyl± 10, taklng/<2 also In the range 
i y c - 10. The results show that Tor a simply supported edge only two or three Iterations 
m ■ 2. or 3) are necessary to obtain solutions of an accuracy of Tour significant digits In 
toe coefficients F n and C n . For a clamped edge one, or at most two steps yield the same 
accuracy, which Is sufficient In most applications. Csrrylng the computation up to m - b In 
a few particular cases shows that an accuracy of six digits can be obtained. We have Indica¬ 


tion. however, that the convergence Is considerably slower Tor values or / In the region of 
the critical load. At the present time, this region la being Investigated with the aim qf 
obtaining approximate values for the buckling load for various values „f the parameter^". 


Hubertus J. Welnltschke 


References: 
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[z] R. M. Simons, "On the Hon-llnear Theory of Thin Spherical Shells," M.I.T. Ph.D. Thesis, 

1955. 


I HOl’ERTIES OF THK SOLUTIONS OF A DIFFERENTIAL EQUATION ARISING IN THE THEORY OF LINEAR 
ELASTIC DEFORMATION 0? SHALLOW SHELLS OP REVOLUTION 


It SB proposed to develop theory and also the necessary tabular Information for use 
In solving for the deflections of shallow shells of revolution. The shells dealt with here 
are parabolic, having as the equation of the undeflected mtdsurface 


H (~) a or H lj S s > 0 


where r Is the radial distance Trom the axis or symmetry and / lr. the height of the mid¬ 
surface above a reference plane. H and a are characteristic lengths. A criterion of shallow¬ 
ness Is that H Is 3mall compared to a. 


Johnson [l] has shown that a solution for the stresses and transverse deflection 
In such a shell satisfies the single (complex) differential equation 




where T la a complex dependent variable who3e real part represents deflection and whose 
Imaginary part represents an Airy 3tress function. Primes Indicate differentiation with re¬ 
spect to r. while the second Independent variable Is the polar coordinate w. The real parameter 
Vc/t> depends on the thickness and elastic constant of the shell. 


Seeking for simplicity solutions for T of the form 


t (r) cos na, 


we obtain an ordinary differential equation for t . 


(&■-■*&- i( i c * zyf-#->.)] -« 


Writing t„ (v) = f m (J> f * ) , f-n H/Z7i> 

and Introducing the new Independent variable 


(M) becomes 


(">) -^ + 2lil=±J * i^s 

+ -ii [>Vn , -W/> f h'(.s-')i\ ^(i) = o 
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Thin equation haa a regular singular point at the origin. Thus at least one poe 
solution aay be found. The lndlclal equation la 

f(x)• (*-»)(*-’) rz 

t [[S-0 (7r-r) - (i i-ZH')] 3*" 


which has the roots 


.. n n 2+n Z-n 

>>■ >• ‘ «’ s’ s 


f(W - (A»-n)O+ n )0»-2-n)O-2+n) 


How, If we write the solution fortsally 




then the A ^ will satisfy a reourslon 


° ' 4fc 

g(» - s z AU-U ♦ s(s-l)A - n 2 (s-l) 


a closed Toro may readily be obtained for A^ In terns of gamma functions: 

( 8 ) . » r (*+<=•- }*) r(Mk- X ) __ 

Here A Q Is not the same as In (7) and ^ 0 , ^ are derined by 

* Jf~ V 's\zs(\-«*i x )+l*»»'- 

I, *■ /jVirO- *h» + 

(We «ay note that. Tor real 8 and real n, >. and X- are also real.) In caBe none of the 
Indices > given In (6) differ by an Integer, four distinct solutions of type (7) will exist, 
one lor each \ values. In the singular cases (which are often of significant interest), the 
method of Probenlus will yield further solutions. 

Much more Information on the characteristics of these functions Is necessary In 
order to solve the boundary value problem. In this project we undertake to find pertinent 
Information such as integral representations, asymptotic behavior, recursion relations, and 
sufficient tabulations to take care of at least simple boundary value problems. 

M. Douglas Mcllroy 

Reference: 

[l] M. Johnson, unpublished paper, M.I.T. Mathematics Department 

ASYMPTOTIC SOLUTION OP A DIFFERENTIAL EQUATION 

The last quarter was spent In research on a problem proposed by Professor C. C. 

Lin, which concerns itself with the study of certsln special functions required for the 
asymptotic solution of differential equations with "turning polnta." In particular, the func¬ 
tions defined by the equations 

(0 ♦ > 2 (*u # + pu 7 ♦ qu) - 0 p,q constants (2) u 7// ♦ 2u 1 ♦ pu ■ 0 

are now being Investigated. 

Tabulation of these functions Is necessary for uniformly valid differential equations 
involving turning points. These were chosen for tabulation since they are the simplest of the 
type desired and can be treated by the method or Laplace transformation to yield explicit 

solutions. 


Joseph Hershenov 
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EVALUATION OP MULTI-CENTER INTEGRALS 


The evaluation of multi-center integrals for use In a tight-binding energy band 
calculation or a two-dIntensions 1 graphite crystal model Is now completed. This task has 
proved to be very laborious even with the aid of the Whirlwind Computer. One compensation, 
however, la that the general techniques required are applicable to all the usual molecular 
integrals and use has been made or programs developed by the writer and by several members 
of the Solid-State and Molecular Theory Group. These multi-center techniques nave been 
described by the writer In a rather lengthy report which will be published In the Quarterly 
Progreso Report of the Solid-State and Molecular Theory Group, April 15 , 195 b. 

Fernando J. Corbatd^ 


COULOMB WAVE FUNCTIONS 


The Integral representation method of computing the Irregular Coulomb Wave Function 
g.(^,>)» has given satisfactory results for the region of Interest indicated :n Fig. 1. The 
particular representation discussed in the lsst few progress reports has been discarded in 
favor of 

yi>»- /vl 

do ' * J l 

O 9 

which does not display the infinite frequency oscillation asa> elated with the former near 
J - 1. The above representation had been tried previously using a rather clumsy nine-point 
numerical Integration technique. We are now using the more efficient Gaussian quadrature 
method. 

To Insure roughly five-place accuracy, the first integral is cut off at f *2 ft* 

and the second Is Integrated up to J - * after which the Integrand Is replaced by tne 
approximate form 

V. {o fTTif - 2>?} 

which can be Integrated exactly from J * 4 to f * <*» . Sixteen point quadratures are used 
over regions ranging from 0.5 to 1,0 depending on how fast the argument of the sine function 
Is changing. 

It is very easy to generate many g »s almultane usly. tf the spacing In tne tables 
Is to be Aand A^» , we calculate for each value off involved in the quadrature, the func- 

c- (Ar>s ; 

It Is then a trivial matter to compute the integrands for the desired mesh of volu«?s. 

In this way the time needed t> compute a g and its derivative may be reduced to about twi 
seconds. If a separate Integration were needed for each point t< be tabulated, the time re¬ 
quired would be about 25 seconds. 

It is hoped that the full production of the Irregular function will start about the 
middle of April and be completed by the end of May. 



I A. / 6 t* a /*, !( If a- * ■** U *+ /• * l ** H 


Arnold Tub la 
Aaron Temkln 
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PIRST APPROXIMATION SOLUTION ON ORE BODY 

During the past quarter aufriclent progress has been made In the flrat mpproxlms- 
tlon solution to merit the continuation of this method of attack on the more general problem 
of the interpretation of field resistivity measurements. 

As outlined in the previous report, a solution for the two-dlmenalonal case mas to 
be developed and a comparison was to be made with the results of resistivity modeling. This 
work is being done In the Department of Geology and Oeophyslcs In conjunction with the Ph.D. 
thesis of Nr. P. G. Hallof. A program for the theoretical solution was written and run Tor 
variations In the depth, shape and resistivity contrast of the ore body which corresponded 
to those cases already modeled. The theoretical results showed a general difference In 
character Tor a vertical and a horizontal tabular ore body, but the comparison with the modeled 
results was somewhat poor. The character was such that the theoretical results from a horizon¬ 
tal ore body corresponded better with those from the vertical ore body In modeling and vice 
versa. Numerical agreement was not satisfactory, and further consideration or the two- 
dimensional case ha9 ceased. 

The investigation was then directed towards the case of three-dimensions. In which 
modeling has also been done by Nr. Hallof. The analysis of this problem la similar in approach 
to that of the two-dimensional case and will not be given here [l] . As a simplification of 
the computational procedure, an approximation was made to the actual potential difference 
measured on the surface of the earth by the receiver electrodes. This was to compute the 
rate of change of potential with respect to the distance from the sender electrode at a point 
midway between the receiving positions and then multiply by the spread separation distance in 
order to obtain the measured potential dlatance. A program to compute the potential differ¬ 
ence as specified above was written and results obtained for variations in the shape, depth 
and contrast of the ore body. The results of one of the theoretical solutions is plotted In 
Figure B, the actual value plotted being the percentage difference In the potential measured 
as caused by the presence of the ore body. The results are plotted In the manner described 
in the previous report, for the case of a horizontal tabular ore body one unit deep, with 
the dimensions of unit thickness, four units wide and eight units long. The location of the 
ore body is shown dotted, and the contouring la done at levels of 0, 10, 20 and AO per cent. 

The results from the modeled case are plotted In Figure A, and it Is seen that the 
comparison between the results plotted In Figure B are quite good. The character of the two 
Is almost Identical, and the numerical agreement is relatively good. However, there are a 
few critical positions of the sender and receivers that give poor results and further work la 
planned to Improve these results. It was expected that some Improvement might be possible 
by the use of discrete receivers In the theoretical solution, rather than the values obtained 
from the derivative of the potential. Another program to calculate this was developed and 
the results are plotted In Figure C. The comparison with the modeled values is somewhat 
poorer In general, but It has not yet been definitely established as to which receiving method 
Is best In general. 

The final goal of this Immediate Investigation Is to determine the optimum alze 
and shape of the ore body, that will give the best comparison with the modeled results. 

With this unit as a building block, It lo proposed to consider any large ore deposit as 
being composed of a number of these basic elements, each with different conductivities. 

The conductivities determined would be those that would best fit the field measurements. An 
economic evaluation of the prospect could then be made by the geophysicist with a knowledge 
of the conductivities, somewhat dependent on the amount of metal present, and of the geology 
of the area. 


Norman P. Ness 

Reference: 

[l] Homan P. Neas, Machine Methods ofComputatlon and Numerical Analysis, Quarterly Progress 
Report Ho. IB, December 15 (1955) p. 21 . * * 
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FIRST APPROXIMATION SOLUTION ON ORE BODY 

During the past quarter sufficient progress has been made In the first approxima¬ 
tion solution to merit the continuation of this method of attack on the more general problem 
or the interpretation or field resistivity measurements. 

As outlined In the previous report, a solution for the two-dimensional case was to 
br developed and a comparison was to be made with the results of resistivity modeling. This 
w.rk Is being done in the Department of Geology and Oeophyslca In conjunction with the Ph.D. 
thesis of Nr. P. Q. Hallof. A program for the theoretical solution was written and run for 
variations In the depth, shape and resistivity contrast of the ore body which corresponded 
to those cases already modeled. The theoretical results showed a general difference In 
character Tor a vertical and a horizontal tabular ore body, but the comparison with the modeled 
results was somewhat poor. The character was such that the theoretical results from a horizon¬ 
tal ore body corresponded better with those from the vertical ore body In modeling and vice 
versa. Numerical agreement was not satisfactory, and further consideration of the two- 
dimensional case has ceased. 

The investigation was then directed towards the case of three-dlmenslona, in which 
modeling haa also been done by Mr. Hallof. The analysis of this problem Is similar In approach 
to that of the two-dimensional case and will not be given here [l] . As a simplification of 
the computational procedure, an approximation was made to the actual potential difference 
measured on the surface of the earth by the receiver electrodes. This was to compute the 
rate of change of potential with respect to the distance from the sender electrode at a point 
mldwey between the receiving positions and then multiply by the spread separation distance In 
rder to obtain the measured potential distance. A program to compute the potential differ¬ 
ence aa specified above was written and results obtained for variations In the shupe, depth 
and contrast of the ore body. The results of one of the theoretical solutions Is plotted In 
Figure B, the actual value plotted being the percentage difference In the potential measured 
os caused by the presence of the ore body. The results are plotted In the manner described 
In the previous report, r t the case of a horizontal tabular ore body one unit deep, with 
the dimensions of unit thickness, four units wide and eight units long. The location of the 
ore body Is shown dotted, and the contouring Is done at levels of 0, 10, 20 and A0 per cent. 

The results from the modeled case are plotted In Figure A, and It Is seen that the 
comparison between the results plotted In Figure B are quite good. The character of the two 
Is almost Identical, and the numerical agreement Is relatively good. However, there are a 
few critical positions of the sender and receivers that give poor results and further work is 
planned to Improve these results. It was expected that some Improvement might be possible 
by the use of discrete receivers In the theoretical solution, rather than the values obtained 
from the derivative i the potential. Another program to calculate this was developed and 
the results are plotted In Figure C. The comparison with the modeled values Is somewhat 
poorer in general, but It has not yet been definitely established as to which receiving method 
Is best In general. 

The final goal of this immediate investigation Is to determine the optimum size 
and shape of trie ore body, that will give the best comparison with the modeled results, 
with this unit as a building block, it 13 proposed to consider any large ore deposit as 
being composed of a number of these basic elements, each with different conductivities. 

The conductivities determined would be those that would best fit the field measurements. An 
economic evaluation 1 the prospect could then be made by the geophysicist with a knowledge 
! the conductivities, somewhat dependent on the amount of metal present, and of the geology 
of the area. 


Norman F. Mess 


Reference: 


[l] Honran P. Ness. Machine Methods ofComputation and Numerical Analysis, Quarterly Progress 
Report No. iti, December 15 (1955) P. Zl. 
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AN APPLICATION OP MONTE CARLO METHODS TO NEUTRON DIFFUSION 


Introduction 

In the last report several methods were described to obtsln approximate answers In 
neutron diffusion problems. The reasons for the selection ol a Monte Carlo method were given. 
This report will give an overall description of the application of this method and an outline 
of the logical design for a program In which this method la used. 

Outline of the Method 

The main feature or the Monte Carlo method is that the behavior of Individual 
neutrons is studied. Thin behavior Is then averaged In such a way that the desired results 
can be obtained; e.g., if we want to know the neutron balance or a reactor, we record the 
number of neutrons that la absorbed or leaked out, and the number of neutrons produced by 
fission. More Information, like the coordinates or each absorption or the neutron fluxes In 
the reactor, can be obtained if that Is desirable. 

The behavior of the Individual neutrons is studied In this way: We derine a ’’bank" 
which la a part of the memory of the computer In which a Bupply of neutronB with their space 
coordinates Is available. If new neutrons are formed by fission, they are put In this bank, 

>r If neutrons are needed, they are taken from thlB bank. 

We start the problem by taking a neutron from the bank. The following description 

can be followed on the attached flowsheet. First we determine now In which direction the 

neutron will fly by calculating from random numbers a set of three direction cosines Buch 
that the direction of the resulting vector Is random In space. As neutrons are produced at 
high energies, we set the velocity switch on fast. The use of this switch Is explained later 
on. We choose now a velocity from the fission spectrum and determine In which velocity-group 
this velocity Is. A velocity-group Is characterized by a number which the computer uses to 
determine which constants to use. E2ch velocity Interval between two chosen limits has Its 
own velocity-group number. 

Now we determine In which part of the system our neutron Is. With this knowledge 
«nd the velocity group, we determine which nuclear constants to use. On the basis of these 
constants and a random number, we calculate how far the neutron will fly. We can now deter¬ 
mine If the neutron at the end of Its flight was still In the same material. If this Is not 

the case, we set the domain switch on out , advance the neutron to the boundary and proceed to 
find out the nuclear properties on the other side of the boundary. If the neutron was still 
In the same material, set the domain switch on Jn. We let a random number decide If the event 
at the end of the path was a capture, flBSlon or scattering. For simplicity only these three 
processes are mentioned. 

If the neutron was captured, we record Its death coordinates for later print out 
and take a new neutron from the bank. If the neutron caused a fission, two or three neutrons 
(av. 2.*»0) will be put in the bank with the coordinates where the fission took place. We 
then take a new neutron from the bank. If the neutron Is scattered, we determine by which 
it'im It was scattered. With this knowledge we are able to calculate the new direction cosines 
i* the neutron by assuming random scattering in the center of mass system. 

If the velocity switch Is on fast, we calculate the new velocity, determine the 
velocity group, and store the constants that belong to this velocity group. 

If the new velocity group Is the thermal group, set the velocity switch on thermal. 
This has the eifect that from that moment on no velocities are calculated anymore, this saves 
time. 


If the velocity switch was on thermal, we do not calculate the new velocity because 
we assume that all thermal neutrons have the same average speed. In this case we calculate 
directly the distance of flight without calculating the velocity. 


we have now completed the cycle. The 
other until It is finally absorbed or leaked out 


me neutron win nave 


«•, v ni T ^ C dar J* arro " 8 * n alowsheet Indicate the part of the program that Is used after 
® c 1 collision. The semi-dark arrows show the part that Is only used several times In the lire 
i i»S#- n ! Utr0n: 4 l i« iS U8ed ° nly after boundaries are crossed or velocity groups are changed. Th' 
light arrows Indicate the part or the program that is used only once In the life of a neutron, 

after its absorption or leakage. 






Conclusions 

All routines outlined In the previous pages hve been programmed, coded and tested 
.•parately. During the next quarter these routines will be assembled and the complete program 
will be tested. Then a number of problems or which the analytical solution Is known will be 
run. Only after all these tests have been completed successfully, actual problems will be run 

Marius Trooat 
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RESPONSE OP A SINGLE STORY REINFORCED CONCRETE BUI1DING TO DYNAMIC BLAST LOADING 

The problem Is to study the dynamic response or the roof system of this building. 
This is a system of two degrees of freedom consisting of a slab supported on relatively 3tlft 
beams. The acceleration of the beam has a great effect on the slab acceleration. 

The loading portion of the problem also has distinctive characteristics since the 
.1 lab area 13 small relative to total roor area and greatly effected by edge vortex action. 

The program to aolve this problem has three main part3: 

(1) The loading portion of the program Is written and is about to be "debugged." 

(2) The integration portion uses the backward dtrrerence technique. It Includes 
calculation of resistances and has been written and "debugged." 


(>) The control portion la almost completed. 


Ralph G. Gray 


RESPONSE OF A FIVE-STORY BUILDING TO DYNAMIC LOADING 

This Is a program to determine the response of actual multi-story rectangular 
rigid-frame buildings aubjected to blast loads. The program as a whole resolves Itself Into 
four parts: 

(1) Control Program 

This program, currently being written, varies the fundamental parameters of 

load. 

(2) Loa d Pro gram 

Written and "debugged," this program computes tne loads due to blast pressures 
on the inner and outer surfaces of the building as functions of time. 

(^) Integra tIon Program 

Thin p- gram has been written. It finds the lateral displacements of the 
floors by numerical Integration, using tne "backward difference" technique to 
determine peek displacements and compares them with allowable values. 

(A) Res lstance Program 

Alan written, this program utilizes the displacements determined in the Integra¬ 
tion program to compute resistances within each story. 

This last program solves a matrix composed of the slope-deflection equations to 
determine Joint rotations and uses the rotations at each end of a particular member to ob¬ 
tain a value cf elastic moment at each end. This moment Is compared with the value at which 

trie member enler3 tne plastic range. If It haa entered this range, the matrix Is modified 
icoord 1ngly. A member changing from plastic to elastic also causes proper modification of 
the matrix. The sequence elaatlc-plastlc-elastlc may be repeated any number of times, values 
of residual plastic strain ("kinks") being accumulated during the process. 

Prom the moments exist rig at the Joints, the story shears are computed and treated 

by the Integration program on the next cycle as resistances. 


Ralph 0. Gray 


GROWTH OF FATIGUE CRACK.'! 


Recent work Is being directed toward investigating the growth of a crack In an 
axially cracked bar loaded In torsion. This work 1 b confined to consideration of bars 
wh ch are perfectly plastic and which are strained to an extent that all stresses are in 
the plant c region. The assumed mode of crack growth Is as follows: The bar Is twisted 
unt. 1 the atraln at a distance In front of the crack reaches the value the crack 
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then advances a distance d f, relieving the stresses; the .ingle of twist, », is increased un¬ 
til the strain x from the crack tip reaches / and again the crack advances, etc. By 
combining this criterion for growth with the equation expressing the strain distribution 
head of the crack l l] , one may express the Inverse of the non-dlmens luna 1 Ued rate r crack 
growth In a rectangular bar as the following Vulterru Integral equation. 

(WL ■ 4r[' - iftl * ft /V ftel W'i 

o 

Numerical solutions to this equation were carried out manually, using two methods: 
The integral was approximated by the trapezotdal rule and a stepwise solutl >n htalned (2j. 
This s ilutlon was used as the first approximation fur the second method--an iterative solution 
. . ng Simpson's rule to approximate the Integral [2]. One Iteration sufficed to show that 
the two solutions were in sufficiently clone agreement. The equation has been programmed 
for the Whirlwind CS computer, using the first method outlined above. 


Joseph B. WaIsn 


References: 


[l J H. A. McCllntock, "Growth of Fatigue Cracks under Plastic Torsion," unpublished paper 
.'submitted to International Conference on Fatigue ol Metals. Institution of Mechanical 
Engineers, London, (195L). pp. 1-12. 

Fzl L. Collatz, Numerlsche Behandlung von Different la lgle Ichungen. Springer, Berlin, (19*i>l). 
1 pp. A11-U15. - 

CONDENSATION OF VAPOR IN A VERTICAL TUBE 

The equations given In tne past report for tne film condensation of a vapor In a 
vertical tube In the ca3e when the vapor la In turbulent flow and tne condensate film In 
laminar flow have been changed. It was found that the ordinary friction factor equation for 
turbulent, single phase flow In a tube could not be used for calculating the anear force Im¬ 
posed by the vapor on the condensate at the vapor-condensate Interface. Data lor the fric¬ 
tion factor In a condenser tube have been published by 0. r. Bergelln, . K. Kegel, F. G. 
Carpenter, and C. Gazley, Jr., In the Journal of the Heat Transfer ana Fluid Mechanics In¬ 
stitute, 1949. These data Indicate that the friction factor of the gas In the condenser 
tube can be several times the value as given by the ordinary single pha3e equation at the 
corresponding gas Reynolds number. The friction factor In the above publication has been 
plotted as a function of the gas Reynolds number and a parameter, Sw-P/sp. 

Using these results, the form of the equations governing the laminar condensate- 
turbulent vapor case did not cnange* although It was necessary to modify the definition of 
the dimensionless variables and the parameters, and to Introduce one new parameter. 

The turbuLent condensete-turbulent vapor case to being analyzed In a manner very 
similar to the laminar condensate - turbulent vapor case. The Nlkuradse velocity dlstrtbut Ion 
has been assumed for the condensate film near the wall while there is still some question 
whether this distribution will be valid near the vapor-condensate interface. The Martinelil 
result to the heat transfer-momentum transfer analogy 13 being used for determining the 
resistance to heat transfer of the turbulent condensate film. 

The laminar condensate-turbulent vapor equations are being programmed for the 

computer. 

Jukka A.Lehtlnen 
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Part 11 

Project Whirlwind 
1. REVIEW AND PROBLEM INDEX 


This report covers the specific period of December 14, 1955 to March 20, 1956. 
During this time, 87 problems made use of 405.8 hours of the 545.6 hours of Whirlwind 
computer time allocated to the Scientific and Engineering Computation (SfcEC) Oroup. 

The remaining 139.8 hours of the allocated time were used for terminal equipment test¬ 
ing and calibration, demonstrations, and various inter-run operations not logged to 
specific problems. 

The 87 problems run during this quarter cover some 16 fields of applications. 
The results of 24 of the problems have been or will be Included In academic theses. In 
these 24 problems, there are represented 23 Doctoral theses, 3 Naval Engineer's, 1 Civil 
Engineer's, and 4 Master's. Thirty-five of the problems have originated from research 
projects sponsored at MIT by the Office of Naval Research. 

Two tables are provided as an Index to the problems for which progress reports 
have been submitted. In the first table the problems are arranged according to the 
field of application, and the source and amount of time used on WWI Is given. In Tab’e 
2-11, the problems are listed according to the principal mathematical problem involved 
In each. In each table, the letter after the problem number Indicates whether the 
problem Is for academic credit and whether It Is sponsored. The code is explained on 
page 28. 
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WHIRLWIND CODING AND APPLICATIONS 


2. WHIRLWIND CODING AND APPLICATIONS 


2.1 Introduction 

Progress reports as submitted by the various programmers are presented in 
numerical order In Section 2.2. Since this summary report presents the combined efforts 
of DIC Projects 6345 and 6915, reports on problems undertaken by members of the Machine 
Methods of Computation (MMC) Group have been omitted from Section 2.2 of Part II to avoid 
duplication of Part I. For reference purposes, a list of the MMC Group problems appears 
below. 


Letters have been added to the problea numbers to Indicate whether the problem 
is for academic credit and whether it ie sponsored. The letters have the following 
significance. 

A implies the problea is NOT for academic credit, is UNsponsored. 

B implies the problea IS for academic credit, is UNsponsored. 

C implies the problea is NOT for academic credit, IS sponsored. 

D implies the problem IS for academic credit, IS sponsored. 

N Implies the problem is sponsored by the Office of Naval Research. 

L implies the problea is sponsored by Lincoln Laboratory. 

The absence of a letter indicates that the problem originated within the SAEC Group. 


List of Machine Methods of Computation Group Reports 

The following problems used Vhirlwlnd I computer time during this quarter, 
but are reported on in Section 2.2 of Part I. 


122 N. COULOMB WAVE FUNCTIONS 

203 B.N. RESPONSE OF A MULT I-STORY FRAME BUILDING UNDER DYNAMIC 
LOADING 

231 B,N. REACTOR RUNAWAY PREVENTION 

329 N. FIRST APPROXIMATION SOLUTION ON ORB BODY 

337 N. NON-LINEAR SECOND ORDER DIFFERENTIAL EQUATIONS IN THE 
THEORY OF ELASTIC SHELLS 


A. Temkln 
A. Tubls 

R. G. Gray 

M. Troos t 

N. F. Ness 

H. Welnltschke 


100 COMPREHENSIVE SYSTEM OF SERVICE ROUTINES 

Upon the dissemination of DCL Memo 117 on a proposed post-mortem program which 
would print out Instructions only, but which would glvw the reference addresses relative 
to floating addresaea and place flad taga In the print-out, the ataff decided that the 
preaent Generalized Poat-Mortem program ahould be revlaed to Incorporate theae feature! 
and to add eeveral others to make the GPM more useful as a means of getting output. The 
goala for the GPM are to give prlnt-outs of storage that are as close as practicable to 
the Input language to facilitate checking and prlnt-outs of data In an easy to use form 
and In an automatic operation that will not uae up very much time. 

A new GPM has been coded but has not yet begun to be checked out. It has these 
facilities: 

(1) The programmer may specify up to 7 different (i)8TOP Instructions as Inter¬ 
mediate break points. These STOP Instructions will not stop the computer or cause any 
logging to take place, but will transfer control to the GPM which will execute a special 
request pertaining to the particular break point and return control to the next Instruc¬ 
tion after the break point. A final break point 00 should be specified. When It Is 
reached another post-mortem will be executed, the STOP will be logged and the computer 
stopped (or the next director tape word read). If the computer stops In any other way, 
the operator nay assume that there was trouble and he will be able to call In the GPM 

to execute a special "in case of trouble" post-mortem by setting the right MIR to 1.00000 
and pressing the Examine and Read-In buttons. The operator may call for any break point 
post-mortem by setting the right MIR to 0.00033, the left M1H to O.OOOOn, where n Is the 
break point number, and pushing both activate buttons and the Read-In button. (Unless 
the Erase button is pushed, the left MIR will be considered to be zero.) The operator 
may obtain scope post-mortems and programmed arithmetic post-mortems (PAPMs) as at present 
(rm4, rl will get a PAPM on the scope), but for more automatic operation It will be 
possible to request these on an fp tape. 

(2) In general, the new prlnt-outs will have the same form as the old. However, 
STOP Is given for sp25 as a Whirlwind instruction and the interpreted instructions 1STOP, 
IN, and OUT are printed. In the flad modes automatic in-out requests will appear as 
AIORQ 0 as explained in Memo OCL-117. Whenever flads are requested, flad tags are placed 

In the body of the prlnt-out. The absolute addresses are always printed in the left-most 
column, but these will usually be the current addresses (the place In core memory these 
words would occupy) rather than drum addresses (the actual locations on the drum). The 
drum address of the first word of the block Is always given on the first (subtitle) line, 
even when data is printed in the format mode where no current addresses are given. There 
will be 16 modes: gda, gdf, of a, off, dla, dlf, dfa, dff, wwa, wwf, wla, • it, iia, Ilf, 
iwa, and lwf , where the last letter calls for absolute or floating address and the first 
two letters call for numbers In the same way as the old GPM and for Instructions the 
first letter sets the mode and the second sets the mode after IN or OUT, l.e., wl and iw 
imply switching.* For numbers, a special format word In place of the Initial current 
address will set the number per line and number per block. Except In this mode, large 
blocks of zeros will be suppressed. 

(3) The punch mode may be selected by an fp tape. The new punch mode will 
give an fb tape complete with title and sp block, but only the 100-octal block length 
will be retained. 

(4) Because flads will be allowed In the addresses appearing In the fp tape 
and the flad table must be preserved for the post-mortem In many cases, the fp tape must 
be read in Immediately after the fc tape containing those flads. As the read-in of tape 
does not result in Immediate execution of the request, but only In storing the request 
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on the butler dru», It will be possible to check for many error, and reject Illegal tape, 
before the prograa 1. operated. 

(S) A PAPM will not be given unlee. It 1. apeclflcally reque.ted at a break 
point. It will be given automatically In ca.e of trouble. 

It 1. expected that the new GPN will be operating routinely In the next quarter. 

A new aeno and/or supplement to the CS Manual will be published to acquaint tho.e concerned 
with the new GPM. 

M. R. Weinstein 

Digital Computer Laboratory 

• ummo DCL-117 FI ad Post-Mortems and ff Tapes, M.R.Weinstein, Hecwber 19, 1955 


101 N. OPTICAL PROPERTIES OF THIN METAL FILMS 

One successful run on sone samples vastly different from ones previously used 
was preceded by a few unsuccessful ones to bring performance and tape requests up to 
date after sore than a year's Interval. 

A. L. Loeb 

Lincoln Laboratory 

120 B,N. THE ABROTHKRMOPKKSSOR 

The present program was found to be unsatisfactory for certain Initial con¬ 
ditions because some Iterative procedures did not converge. Small changes should once 
more make the present program satisfactory. 

Future calculations are planned to check the experimental results that are just 
being obtained in the Gas l\irblne Laboratory. 

A. J. Erickson 

Mechanical Engineering Department 

126 C. DATA RKDUCT10N 

Problem 126 Is a very large data-reduction program for use In the Servomechanisms 
Laboratory. The overall problem is composed of many component sections which have been 
developed separately and are now being combined into complete prototype programs. Desc¬ 
riptions of the various component sections have appeared in past quarterly reports. 

After the development and testing of the prototype Whirlwind programs Is completed, the 
programs will be re-coded for other, commericaily available, large scale computers, 
(probably the Univac Scientific 1103, IBM 7Q1 and IBM 704 computers), for use by Interested 
agencies for actual data reduction at other locations. The programs are currently being 
developed by Douglas T. Ross, David F. McAvinn, Walter K. Weiaablum, Benson H. Scheff, and 
Dorothy A. Hamilton, Servomechanisms Laboratory staff members with the assistance of 
John F. Walsh. This work Is sponsored by the Air Force Weapons Guidance Laboratory. 


reason extensive use is made of output oscilloscopes so that the computer can communicate 
with the human, and manual Intervention registers so that the human can communicate with 
the computer In terms of broad ideas, whlls tbs computer is running, and hav# the com¬ 
puter program translate these ideas into ths detailed steps necessary for program modi¬ 
fication to conform to ths human operator's decision. The program which does this trans¬ 
lation and modification Is called the Manual Intervention Program (M1V). The most recent 
version of ths prototype data reduction program Is called the Basic Evaluation Program. 

The basic keatures of the M1V Program were described in the previous progress 
report. Uirlng the past quarter the Basic Evaluation Program was completely debugged 
using the Mistake Diagnosis Routine (described earlier) to extract pertinent information 
at over seventy break points. The nsture of the Evaluation Program Is such that this 
debugging process would have been virtually impossible without the use of the MDK program. 
In addition, an entire new section of the Evaluation Program, which will allow automatic 
starting of the data reduction process at any point, was incorporated and checked. 

The decision to operate all Problem 126 programs under Director Tape Control 
has materially Increased the effectiveness and reliability of these programs. The re¬ 
coding of these programs for the Univac Scientific 1103 Computer using the facilities 
of the WWI-1103 Input Translation Program of Problem 256 has begun in earnest. 

The Logging Program section of the M1V Program is now undergoing changes to 
provide greater flexibility of operation. The function of the Logging Program is to 
record logical Information on a magnetic tape unit, giving all pertinent information 
concerning all manual interventions, l.e., which buttons have been pushed and what cer¬ 
tain parameter values are, and similar information concerning all displays seen by the 
human operator. The complexity of the manual operations makes such a complete record 
essential to the successful use of the programs. 

The logical magnetic tape output of the Logging Program is to be processed at 
the end of a run by the Editing Program which will give a delayed printer tabulation of 
the log In English phrases and numbers of appropriate types. A preliminary version with 
very restricted capabilities is now in use, and the elaborate complete version is now 
being written. Mr. Frank C. Helwlg of the MIT Digital Computer Laboratory is a major 
contributor to the Logging and Editing Programs. 

A very flexible and elaborate Scope Plot Program is also being written. This 
program will give axes calibrated with four sizes of pips arranged in a pattern similar 
to that on a ruler, a faint grid In the background, plotted points connected by faint 
straight lines, alphanumeric labeling and automatic Initiation of new frames whenever 
any boundaries are exceeded. The program is fed a series of parameters which tell the 
desired format, scaling, and labeling and then is fully automatic. The program can be 
initiated either automatically or through the MIV program. 

Both the Editing Program and the Scope Plot Program require the use of arbit¬ 
rary sequences of alphanumeric characters. For use in these and future programs a pair 
of general programs are being written. The first program accepts arbitrary six-bit 
characters from Flexo Tape and alphanumeric code words in a special language, and stores 
them in arbitrary locations with optimum packing. The second program unpacks the code 
words and characters, and under their direction distributes the arbitrary Information to 
selected in-out units. 


The nature of the problem requires extreme automatlclty and efficiency In the The Unpacking Program is being realized by programming a new type of computer 

actual running of the program, but also requires the presence of human operators in the which has three modes of operation of its control element. Each instruction in the 

computation loop for the purpose of decision making and program modification. For this program to be executed by this computer lo tagged whether or not It la to be executed in 










APPROVED FOR PUBLIC RE 


WHIRLWIND CODING AND APPLICATIONS 

each of the various modes. A system of Jump and remenber-Jump instructions, which can 
change both location and mode, gives this type of computer extreme flexibility an 
compactness, and should prove to be very educational from all standpoints. The future 
work to be done by this group has a very real need for investigations of this type. 

D. T. Ross 

Servomechanisms Laboratory 

IM S ti £ C SUBROUTINE STUDY 

Subroutines are being checked out for function evaluation, via truncated 
Chebysh«v and Legendre expansions, for root extraction (Barristow and Bernoullli), and 
for matrix converses (Craig’s method and Frame's method). 

M. S. Watkins 

Digital Computer Laboratory 


155 N. SYNOPTIC CLIMATOLOGY 

A complete discussion of the research carried on by the Synoptic Climatology 
Project,under sponsorship of the Office of Naval Research, may be found in the final 
report, "Studies In Synoptic Climatology" published by the MIT Department of Meteorology, 
March 15, 1956. 

Programs for this problem were written by members of the Synoptic Climatology 
Project under the supervision of Dr. T. J. Malone, formerly of the MIT Department of 
Meteorology and Professor E. N. Lorenz, MIT Department of Meteorology. 

This problem is now terminated. 

E. A. Kelley 
Meteorology Department 

1:2 N. ANALYSIS OF A SCATTERING EXPERIMENT IN NUCLEAR PHYSICS 

A fit of various combinations of phase shifts has been obtained at discreet 
points for low energies. Now, a slightly finer mesh of points is being used in the same 
region to determine the set of phase shift combinations which will give the truest physical 
picture. 

Future computations will be based on this set, as it is planned to project out 
computations to try to find a fit at middle and high range energies. 

Programmers working on this problem are E. Campbell and E. Mack of the Nuclear 
Science Laboratory. 

B. Mack 

Nuclear Science Laboratory 

179 C. TRANSIENT TEMPERATURE AND STRESS OF A BOX-TYPE BEAM 

Closer correlation with experimental temperatures was achieved by improving 
the value of contact admittance. Future plans are to improve stress correlation. 

J.C.Lorla 

Aeronautical Engineering Department 
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193 L. EIGENVALUE PROBLEM POR PROPAGATION OP ELECTROMAGNETIC WAVES 

For purpose, of comparison with the results already obtained for a bilinear 
nodel of the atsoaphere at SO Me, we plan to carry out similar calculations for an 
Inverse-square model. This work will Involve computation of Bessel functions of terse 
complex order and argument. 

Calculations of 410 Me and 3000 Me have been continued, as described previously. 
Tapes have been prepared with which eigenvalues, eigenfunctions and mode eusa aay be 
calculated at one time, rather than In two parte, as was done formerly. Preliminary 
tapes have been made for the Inverse-square model. 

The 410 Me and 3000 Me calculations with ths bilinear model will be continued 
until completed. Calculations will be carried out for the Inverse-square nodel at 50 Me, 
for comparison with computations at this frequency which wwrw completed during the 
preceding quarter. 

Programmers working on this problem arm Professor H.B.Dwight end Dr. R.M.Rlng 
of Lincoln Laboratory. 

H. B. Dwight 
Lincoln Laboratory 

194 B,N. AN AUGMENTED PLANE WAVE METHOD AS APPLIED TO SODIUM 

Unfortunately, although the routines mentioned in the last report have been 
Individually tested, the tests of the complete program have not yet been successfully 
completed. 

In the meantime, it has been decided to uee the symaetrlzatlon routines of the 
matrix element generation routine In conjunction with the B ve routine to obtain 

E vs Eq curves for symmetrized combinations of augmented plane waves. 

Detallo can be found In the Quarterly Progress Reports of tbe Solid State and 
Molecular Theory Group, MIT. 

M. M. Saftren 

Solid State and Molecular Theory Group 

216 C. ULTRASONIC DELAY LINES 

This period has been spent In debugging tbe program and In making a few teat 
runs. Tbe program ie very complex, consisting of an Iterative loop which contains a 
linear program. The reason for Iteration with a linear program Is that the restrictions 
(Inequalities) are not linear, but well-behaved analytic functions In many variables, 
for which linear approximations consisting of the first terms ot Taylor series are made. 

The sizes of the problems to be solved vary from 40 restrictions with 20 
variables to 90 restrictions with 32 variables. The tlae needed for one cycle varies 
correspondingly from 3 minutes to IS minutes, most of which Is taken by tbe linear 
program. Usually 3 cycles or less suffice for convergence. 

There are about five of the larger problems yet to be solved and several more 
test runs using the smaller problems will be made to determine the mode of operation 
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preferred. 

The programmers working on this problem are Richard Bishop and John Ackley. 


R. Bishop 

Mathematics Dopartaent 


219 LINEAR PROGRAMMING - TRANSPORTATION PROBLEM 

The Whlrlwlnu program for the solution of the classical transportation problem 
has been completed and proved on a number of trial problems, and is now available for 
general use. The program will handle problems with up to 127 plants and up to a total 
of 401 plants and customers provided the product of the number of plants and the number 
of customers Is less than 7169. Several sets of data for a practical distribution 
problea with 9 plants and 69 customers have been tried. The average computer time for 
each solution was about 1.5 minutes. 


The Classical Transportation Problem 


The transportation problem may be phrased as follows A company operates m 


plants producing a commodity, the 1 


th 


of which can supply units of the commodity. 


The company sells Its production to n customers, the J 


th 


cf which demands D^ units of 


the commodity. The costs of manufacturing and transporting a unit of the commodity from 


plant 1 to customer J Is C 


1J' 


It is desired to find the number of units that should 


be shipped from each plant to each customer, in order that the total cost of the operation 
be a minimum. Thus the problem may be stated mathematically as: 


Minimize 


£ 

i.j 


c ij x ij 


subject to the constraints 

I, . 


^ * 8 i 


x u-° 


( 1 ) 


( 2 ) 


(3) 


(4) 


This is a special case of the general linear programing problea. 
The Stepping Stone Method » 
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to aero in order that equation. 2 and 3 remain satisfied 
.111 again have exactly . + n - 1 non-zero"^ ' 


Thus the nee feasible solution 


reduce C. nl ..t'o^*^ u “^‘aTopti^ "° “ *° 


U 


■ake. it faster than cmpu er r tlchnlquei q for"**!*!° f “° rlng th# fM *“>le solution which 
elsewhere. It 1. planned ^ of o! * Proble. in us. 

Journal of the Association for Cosputlng " t,Chn,q “* b * publlah,d the 

th. MIT « J - B - —»• - ■-». *de. of 


J. B. Dsnnls 

Elwctrlcal Engineering Department 


Dantzlg, 0. B., "Application of the Simplex Method to 
Activity Analyst» of Production and Allocation 
4 Sons, New York, 19S1, pp 359-373.- 


a Transportation Problea", 
T.C.Koopmana, Bd., John Wiley 


Charnes, A., and Coop.r, W. w,, 
Programing Calculations in 
Vol. 1, No. 1, Cctober 1954 


The Stepping Stone Mathod of Explaining Linear 
Transportation Problems", Management Science , 


225 N. NKUTRON-DEUTTiRO 1 SCATTER INC 

The calculation of the neutron-deuteron scattering lengths with th. trial 
function described in previous report, ha. been completed. . u^qu! Iln^ mr. Jound 
or the variational expressions for k cot 8 In both quartet and doublet spin states 

-O detained .how little effect of polarl.atim In^h^m 
intern reran ’ t* * Slt that theae r ® 8ult ® a* - ® Inconclusive since the following two 

r\Tv*" P °” lble <a) the Phy ' lc<11 >• indeed unpolarlx.d 

) the trial function used is not adequate to describe the physical syst.a. 

.. , , In “ affort to thro » farther light on the situation, it 1. planned to continue 

con.ld i ” U * ing ° ther t0rm ‘ tOT th * trlal fl “di°h. Two such form, arc being 
function ^,r e rH° nti i n * tW ° “ <ldlt ‘ c,n * 1 Panameter. to describe th. no-polarlzatlon wav. 
Both of ih ! ? ber contaln ' 4 “ ore flsalbl. description of the polarization 

mlStlL 5«Ha?T W trlal ^ unctlon " ire ai-ilar enough to th. original one so that the 
existing formalism can be taken over with little change. All the required Integrals 

amount Of Ln?" °’' aluated ' * nd the 8xl8tln K Programs can be used with only a small 
amount of modification. This calculation Is now in progress. 


A set of which satisfies equations 2, 3 and 4 is called a feasible solution 

to the problem; if the set also minimizes C, it is an optimum feasible solution containing, 
in general, exactly m + n - 1 non-zero X^. The stepping stone method of solution consists 

(1) Generating a feasible solution having exactly m + n - 1 non-zero x^ . 

(2) rinding a zero X^^ which, if allowed to be positive, would yield a 


L. 8artori 
Physics Department 


• The results of this calculation have been submitted to the MIT Physics Departmint" 
as part of a Ph.D. thesis by the author. 
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22’> D. INVESTIGATION OF TIB VORTIC1TY FIELD IN TIB GENERAL CIRCULATION OF THE ATMOSPHERE 
The following program* hava bean completed during the pant quarter 

a) A 30th order asymmetric matrix Inversion program (24.6 Input). 

b) A matrix transpose program. 

c) A relaxation program for V 0 over a 60 point square grid. 

d) The program for the computation of N(h, A, y) In 

equation (1) j^see Summary Report No. 42^ . 

e) The program for the computation of P . 

Plans for the next quarter Include the writing and checking out of two programs 
for the Relaxation, over two 107-polnt graded net gride, of the operator 

C<y) !t~ * H(y> Jv ( |rT ) ■ 1 ( ^m'^T : y) Tt * 0 

The programmers working on thle problem are Duane Cooley and Martin Jacobs. 

R. L. Pfoffer 
Meteorology Department 

234 N. ATOMIC INTEGRALS 

Subprograms have been written and tested which evaluate the radial Integrals 
arising from operators 

L M i y 

r Y^ and 


U-l 


for all values of L. The functions Y** are spherical harmonics. These programs also 
evaluate radial Integrals for the diamagnetic susceptibility operator 


rV 


and the contact hyperflne Interaction 
6(r) 

2 V 


All of these programs are based on simple recurrence relations for matrix elements of 
these operators between analytic atomic orbitals of the form 


= M 

■ s a 


- Y ar r A+.I 


« Y 0 


where is positive and A Is a non-negative integer. 


R. K. Nesbet 

8olld State and Molecular Theory Qroup 


236 C. TRANSIENT RJISPONSE OF AIRCRAFT TO HEATING 

Ikirlng this quarter the basic temperature distribution computing program has 
been modified In such a way as to permit calculation of the transient temperature and 
stress response of a T section which Is half of a doubly symmetric I section subjected 
to any given continuous timewise variation of the boundary conditions, (l.e., heat 
transfer coefficient h and adiabatic wall temperature T ) 

aw 

Eleven satisfactory runs were made. These runs Included accelerated dives and 
level flight decelerations of Aluminum, Titanium, and Inconel X T-sectiona. The results 
have been Interpreted In terms of structural weight penalties. The results are further 
reported In a paper entitled "Some Structural Penalties Associated with Thermal Flight" 
by Professor James W. Mar and Luclen A. Schmit, which was presented at the Spring Con¬ 
ference 1956 of the Aviation Division of the ASME In Los Angeles, California. 

A study of thermally unsymmetrlcal systems is contemplated for the future. 

Programmers working on this problem are L. A. Schmit and H. Parechanian. 

L. A. Schmit 

Aeronautical Engineering Department 
241 B. TRANSIENTS IN CONTINUOUS DISTILLATION COLUMNS 

A program has been written to determine the importance of the non-linear 
variation of vapor composition with liquid composition leaving the plates in a dis¬ 
tillation column. If a linear variation for each plate Is assumed and the usual 
simplifying assumptions are used, the transient behavior of a column In which a step 
change has been made In one or more of its operating variables can be described by a 
set of simultaneous linear first-order differential equations. The solution of any 
such set of equations can be carried out by standard techniques to find the roots of the 
determinant of the coefficients of the dependent variables and the constants from the 
initial and boundary conditions of the problem. The complexity of the solution Is such 
that a high-speed computer Is necessary to carry out those calculations. 

A comparison of the results obtained in this manner with the results obtained 
earlier by using a finite difference approximation to solve the more accurate non-1inear 
equations should show how Important the non-llnearltles are in many cases. 

S. H. Davis, Jr. 

Chemical Engineering Department 

244 C. DATA REDUCTION FOR X-l FIRE CONTROL 

This problem is concerned with computing from fire control signals how far 
fictitious projectiles miss an observed target. The program tape was ready prior to this 
quarter. fork during this quarter consisted of reducing new target run data by use of 
the program tape. 

Future plans for the problem call for still more data reduction, perhaps also 
using supplementary programs (yet to be prepared) for clearer presentation of the 
reduced data. 

J. M. Stark 

Instrumentation Laboratory 
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245 N. THEORY OF NIUTRON REACTIONS 

Some of the coding has been revoreed to increase the accuracy of the results 
for t ■ 5 and 6. Results have been obtained for 6 * .2 3 * .15 for x's up to 1.6, 

including x * 0. 

Future plans are to revise the progroa so that 5 is a function of x^ Instead 
of a constant. 

E. Campbell 

Nuclear Science Laboratory 
253 N. AN AUGMENTED PLANE WAVE METHOD FOR IKON 

E vs. E curves have been completed for both the face-centered and body- 
centered structures. These curves serve to define the augmented plane waves entering 
Into a variational procedure which in turn gives the so-called K vs. k curves. 

J. Wood 

Solid State and Molecular Theory Group 

256 C. WHIRLWIND I - UNIVAC SCIENTIFIC 1103 TRANSLATION PROGRAM 

The vocabulary of the WW1-1103 input translation program (as described in 
Summary Roport No. 41) is currently being enlarged to include a more generalized form of 
numbers. Programmers will be permitted to write numbers of the form 

1 12.34 • 2 3 • lo” 1 

with or without as many associated factors of powers of 2 and 10 as desired, wherever 
they were previously allowed to write integer numbers. In addition, the post-mortem 
program has beon modified to print an indication of all the illegal words detected by the 
translation program. 

Work la expected to be completed on the translation program during the next 
quarter, and a final report on tho program will be issued as a Digital Computer Laboratory 
memorandum. 

Programmers working on this problem are F. C. Helwig and J. M. Frankovich. 

J. M. Frankovich 

Digital Computer Laboratory 

257 C. HORIZONTAL STABILIZER ANALYSIS 

IXirlng this quarter the program for the solution of the two rigid body equations 
of motion (vertical translation and pitching) has been checked out and is now ready to be 
Incorporated In the overall analy.la of the horizontal atablllzer. The program for the 
regponae of a thirteen lumped maee and a one lumped Inertia system, representing a typical 
present-day fighter aircraft, to a etep-funitlon Input la almoat completely checked out 
In the linear elaetlc regime of deformation. 

In the future, the two programs mentioned above will be Incorporated Into one 
program to yield realistic structural response, or deformations, of a present-day fighter 
to a gust type input. The responae of the fighter will be solved for in tts elastic 
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regime of deformation and also beyond, when the structure Incurs permanent deformations. 
The solution of the differential equations of motion of the system Involves the use of 
numerical Integration In a atep-by-atep procedure. 

8. Crlsclone 

Aeroelastlc and Structures Research Lab 

260 N. ELECTRONIC ENERGY OF THE OH MOLECULE 

Most of the work done In this quarter baa Involved chocking the results 
previously obtained and making use of programs (described In earlier Summary Reports) 
for repeating some of the calculations Involved In calculating the binding energy of the 
molecule at new lnternuclear distances. 

A. Freeman 

Solid State and Molecular Theory Group 

261 C. FOURIER SYNTHESIS FOR CRYSTAL STRUCTURES 

The crystal structure of pectollte mentioned In earlier reports was finished 
during the last month, while wollastonlte Is nearing completion. 

The refinement of diglycine hydrochloride, also explained In earlier reports, 
is finished as far as projections are concerned. Only one more refinement in three 
dimensions remains to be done. 

The Isomorphic diglycine hydrobromide is under Investigation, and refinement 
is carried on In Its main projsctlons by difference syntheses. 

Further progress toward the completion of the structure of Jamesonlte has been 

made. 

Programmers working on this problem are Professor M. J. Buerger, T. Hahn and 
N. Nlizekl of the Geology Department, and J. Roseman of the Digital Computer Laboratory. 

T. Hahn 

Geology Department 

262 N. EVALUATION OF TWO-CENTER INTEGRALS 

The exchange Integrals evaluated by Merryman’s program, was not found to be 
sufficiently accurate at lnternuclear separations of one and two atomic units. These 
integrals are now being evaluated by Corbat^'s program. 

H. A. Agbajanian 

Solid State and Molecular Theory Group 

264 C. OPTIMIZATION OF AIRCRAFT ALTERNATOR REGULATING SYSTEM 

The method for minimizing a function of n constraints In the presence of 
equality constraints by means of the method of steepest descent has beon proven for two 
simple mathematical problems. 

Work is being done on improving the above program, adapting It to large scale 
problems, and modifying It for inequality constraints. 
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Programmers working on this problem are R. R. Brown and J. B. Dennis of the 
MIT Electrical Engineering Department. 

R. R. Brown 

Energy Conversion Laboratory 

266 A. CALCULATIONS FOR THE MIT REACTOR 

During this quarter, work continued on calculations for the MIT Reactor using 
parameters that are closer to the actual reactor design. The same basic two-group, four 
region model Is being used. Other models are being developed and calculations using 
these models will be tried In the future. 

The reactor power transient studies, carried on with Whirlwind I and the 
Dynamic Analysis and Control Laboratory analog computer, have been temporarily suspended. 

M. Troost of the MIT Chemical Engineering Department Is doing the programming 
for this problem. 

T. Cantwell 

Nuclear Engineering Group 

270 B. CRITICAL MASSES IN DgO MODERATED REACTORS 

During this quarter tapes have been prepared so that a large number of cases 
may be run. The following Is a list of all methods now In working order; 

Spherical Geometry 

1. Serber Wilson one-group diffusion. 

2. Three-group Fourier Expansion. 


Cylindrical Geometr 


1. Two-group diffusion. 


2. Three-group Fourier Expansion - both approximate and exact expansion. 

The only remaining program to be debugged Is the spherical multi-group transport. This 
should require only a few weeks. 

Future plans Involve testing the worth of these methods by comparing them with 
each other, and with what experimental data Is available. 

If there Is time, a completely new method Involving expanding the flux In N 
groups by Fourier analysis will be tried. This should be free of certain restrictions 
that hamper the use of the three-group Fourier technique. 


J. R. Powell 
Nuclear Engineering 
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272 L. GENERAL RAYDIST SOLUTION 


A now method of solution has been tried during this quarter. This method Is 
purely algebraic and solves for a set of coefficients for an eight order polynomial from 
the coefficients of the original three second-order non-linear equations. 

This method proves the problem to be well conditioned and shows the prospect of 
solution to the desired accuracy. The program at present is In a debugging stage. 

A. Zabludowsky 

Digital Computer Laboratory 

273 N. ANALYSIS OF AIR SHOWER DATA 

The last progress report describes how the method of steepest descent Is used 
to make a least-squares fit of an empirical function f to air shower data. That report 
mentioned the necessity of speeding up the computation. In order to do this the function 
has been changed from an exponential form to the following; 


This function gives an adequate fit to the data and the time per shower Is only about 
30 seconds. But It has the defect that a good fit usually can be obtained with different 
values of 0. For this reason, we are planning to go back to the exponential function 


The computation time per shower may still be less than one minute since the program Is 
now capable of making good estimates of the initial values of the variables. Also, we 
plan to develop a rapid exponential subroutine. 

The problems connected with the empirical function have not prevented us from 
analyzing some data. So far, we have results on about 50 showers. This work will 
continue along with the above mentioned problem and other testing of the program. 


Programmers working on this problem are Professor G. 
the MIT Physics Department. 


Clark and F. Scherb of 


F. Scherb 
Physics Department 


275 B. BUCKLING OF SHALLOW ELASTIC SHELLS 


The work on the problem of the buckling of a compressed shallow cylindrical 
shell Is virtually complete. The results of this problem and the problem of the buckling 
of a hyperbolic paraboloidal shell loaded by Its own weight will appear shortly In a 
doctoral thesis for the MIT Mathematics Department. 

A. Ralston 

Mathematics Department 
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278 N. ENERGY LEVELS OF DIATONIC MOLECULES (LiH) 

The Ground Electronic State 

The ground electronic state as determined spectroscopically occurs at an 
equilibrium tnternuclear separation of 3.0141 atonic units. The spectroscopic value of 
the vibrational constant, CJ e . ls 1 -‘ 106 * 1()3 c " • and the e*P« rl “ ental valuo ot the 
dissociation energy is 2.5 ± 0.2 electron volts. 

The configuration Interaction problem has taoon carried out by considering only 
Configurations of proper ground state symmetry which occur when the ls shell ls kept filled 
The other two electrons may occupy the lithium 2s and 2p orbitals and the hydrogen ls 
orbital. Six configurations of appropriate symmetry occur. 

The value of the binding energy has been determined for seven tnternuclear 
distances between 2.0 and 6.0 atomic units. The potential energy curve which results has 
Luen approximated by a Morse curve, from which wo obtain a value of 1.212 x 10 3 cm* 1 for 
tie vibrational constant, 3.245 atomic units for the equilibrium tnternuclear distance, 
ad 1.669 electron volts for the binding energy. 

The First Excited Singlet State 

The first excited electronic configuration has the same symmetry as the ground 
stato (*£*). It occurs at an equilibrium separation of 4.9052 atomic units, and lies 
3,287 electron volts above the ground state. The vibrational constant, determined 
spectroscopically, ls 234.4 cm* 1 . 8ince the symmetry ls unchanged in going from the 
ground state to the first excited state, it is reasonable to assume that a determination 
of the next to lowest eigenvalue of the secular determinant might be expected to give a 
fair approximation to the first excited state. 

The calculated potential energy curve ls very broad as expected, with a minimum 
very close to the observed value and about 0.29 electron volts too high. A power series 
form, 

U(r) » »- 0 f 2 il + Sj J + a 2 ? 2 +-..) where f w —-—— 

e 

has been used to represent the potential energy near the minimum. The vibrational 
constant, found from a Q , was in good agreement with the experimental value. 

SCF LCAO-MO Calculation 

We Intend to compare the results of the configuration Interaction calculation 
with a self-consistent field molecular orbital treatment. Adequate programs have been 
developed by Heckler and Nesbet for the Whirlwind computer, and are available for eval¬ 
uating the electronic energy of lithium hydride by this procedure. 

Tills phase of the work on LiH ls being Initiated, and further work ls In 
progress on the electronic binding energy at more extended lntemuclear distances using 
the conventional configuration Interaction scheme. 

A. M. Karo 

Solid State and Molecular Theory Group 


285 N, AUGMENTED PLANE WAVE METHOD AS APPLIED TO CHROMIUM CRYSTAL 

The display routine mentioned In the last report has been tested and will be 
used to examine the discontinuity of slope In augmented plane waves. 

8 ws. E q curves which are being calculated for a slight change In the 

chromium potential indicate the sensitivity of these curves to potential, at least In 
the transition metals, The potential used was obtained In a way not usual for energy 
band calculations. However, other members of the Solid State and Molecular Theory Group 
using Whirlwind I will soon furnish Hartree-Fock wave functions for the transition 
metals. From these It ls hoped that potentials suitable for an energy band calculation 
can be obtained. 

Further calculations await the final testing of routines mentioned under 
Problem 194. 

Details can be found in the (fciarterly Progress Reports of the MIT Solid State 
and Molecular Theory Group. 

M. II Saffren 

Solid State and Molecular Theory Group 

288 N. ATOMIC WAVS FUNCTIONS 

Approximate atonic wave functions, expressed as sums of determlnantal functions 
constructed from orbitals which are linear combinations of basic exponentlsl one-center 
orbitals, have been obtained for He, Li, B, and Be. 

The calculations on Li and B are primarily concerned with evaluating the 
electronic contribution to the hyperfine structure, making: use of auxiliary' programs 
described under problem 234. 

A series of calculations on the isoelectronlc systems Be, LIH, He 2 and KL1 
have been carried out by the one-center method(see Quarterly Progress Report, SSAMT Group, 
April 15, 1956), at a fixed internuclear separation of 3.0 atomic units. The difference 
in total electronic energy of these systems can be obtained from the mean value 



taken about the nucleus with charge Z in this lsoelectronlc system (the charge on the 
other nucleus is 4 - Z) . These calculations are Intended to explore the feasibility of 
calculating molecular binding energies by this Indirect method. 

R. K. Nesbet 

Solid State and Molecular Theory Group 

Approximate atomic wave functions are being obtained for the following transition 

elements: 



These calculations are being done in an attempt to obtain criteria concerning the choice 
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ol ono-center orbitals and the Initial approxlaatlon for tha density natrlx for transition 
elenent calculations. 

R. Watson 

Solid State and Molecular Theory Group 
290 N. POLARIZABILITY KFFBCTS IN ATOMS AND MOLKCULKS 

To describe the distortion of atoas In a uniform field, we have set up the 
perturbed Hartree-Pock equations for a general external perturbing field of arbitrary 
order and symmetry. This atomic polarizability problem corresponds to a first order 
perturbing potential of the form r cos 9. As the first step In solving the self- 
consistent Hartree-Pock problem and also to get functions suitable for describing the 
distortion of an atom In an Ionic molecule we are using an equation derived by Sternhelmer. 
The derivation of the perturbed Hartree-Pock equations, a discussion of the physical 
assumptions and of Sternhelmer's equation has been given in recent issues of the Solid 
State and Molecular Theory Group Quarterly Progress Report. 


where 


Sternhelmer’s equation Is: 



t-wf’ u i*-*.r 
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* The desired perturbed function characterized by ft and J? '. 

U * The known unperturbed atomic function characterized by 9 and given as a sum 

of terms each of which Is of the form r n e~ ar . 




With A. Temkin (Morse Group) we have set up the solution of this equation on Whirlwind. 

At present we have obtained solutions for the I ' « 2, Jt = 1 case for the F", Nj , A , 
and Nj ions (the various ions being characterized by the different U 's) by an Inward 
integration technique using the Whirlwind Library Subroutine DE 2 Kutta-Glll. 

Besides this case there are also important solutions corresponding to Q ' * 0, 

C * 1 and 9 ' » 1, 9*0. For these we have used an outward Integration also with the 
Kutta-Glll subroutine but so far we have not been able to get satisfactory results either 
with respect to the form or magnitude of the functions. In particular, the number of 
modes for the case 9 ' “0,9 * 1 has changed as we have changed our aeah size. 
Unfortunately, in contrast to the case of the homogenous ordinary differential equation, 
there appears to be no general theory 2 relating to our unhomogenous equation. The 
equation possesses great sensitivity to the mesh size since we have a very rapidly vary¬ 
ing one at large r. This has necessitated several revisions in both the number and 
distribution of our mesh points and at present we are investigating a procedure employing 
160 points and 19 different mesh intervals. 


L. C. Allen 

Solid State and Molecular Theory Group 
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293 C. ROLLING BEARINGS 

Paat reaulta wore unaatlafactory. Some reaearch la preaently being done on 
a new approach to the solution of the non-linear simultaneous equations. 

Future attempts may be made If a sufficiently promising method la davsloped. 

A. Shashaty 
Lubrication Laboratory 


297 B. DIFFUSION BOUNDARY LAYER 

A general description of the problem was given in Summary Report No. 43. The 
computations completed during the past quarter were concerned with a revised form of the 
energy equation. Previous results from iteration programs for the concentration and 
momentum equations were employed to obtain the additional solutions. The program will 
be concluded during the next quarter. 


J. Baron 

Naval Supersonlcs Laboratory 


300 L. TROPOSPHERIC PROPAGATION 

Values of the Fresnel Integrals for complex arguments have been computed using 
convergent power series. Satisfactory values have been obtained for constant imaginary 
part of the argument as a function of the real part over ranges in which the convergence 
is not too slow. 


J. F. Roche 
Lincoln Laboratory 


306 D. SPECTRAL ANALYSIS OF ATMOSPHERIC DATA 

During the past quarter the writing of the program was successfully completed. 
It is expected that, in the near future, extensive computations of the type outlined in 
Summary Report No. 43 will be made using this program. 

B. S&ltzman 
Meteorology Department 

309 B,N. PURE AND IMPURE POTASSIUM CHLORIDE CRYSTAL 

The calculation of molecular KCl Integrals outlined in previous reports has 
been completed, and the combination of these Integrals into matrix elements Is nearing 
completion. Since the last report, however, It has been found that second neighbor 
Cl-Cl matrix elements must also be Included. This correction Is now being made. It is 
expected that the main objectives of the calculation will be attained during the next 
quarter. 

L. P. Howland 

Solid State and Molecular Theory Group 
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312 L. ERROR ANALYSIS 

This problem, which was described In Summary Report No. 43, has been In 
working order for several months now, and Is being used frequently for production runs. 
Modifications are being prepared to make the program represent more accurately the real 
physical situation. 

L. Peterson 

E. Hutcheson 

Lincoln Laboratory 

314 C. HIGH ORDER POLYNOMIAL FACTORIZATION 

The present program has been extended to factor a series of polynomials on one 
tape. This resulted In a substantial saving of total problem time. Thus, the tapes 
fc 141-227-79 and fc 141-227-1 were replaced by the binary tape fb 314-88-4. A typical 
result is: factoring ten 10 th order polynomials In 2.5 minutes. This permits 

economical use of Whirlwind In computing root loci for systems stability analysis. 

Y. Shulman of the Aeroelastlc and Structures Research Laboratory has been the 
programmer for this section of the problem. 

V. w. Howard 

Y. Shulman 

Aeroelastlc and Structures Research Lab 

317 C. EXTRACTION OF STABILITY DERIVATIVES FROM FLIGHT TEST DATA 

The aerodynamic terms of the equations of motion of a high-speed aircraft have 
been calculated using the method reported In NACA TN3288 (December 1954). It has been 
determined that the type of simulated response used In the study to date adversely 
affected the accuracy of the results. A program has been written, therefore, which will 
provide simulated responses to controlled inputs. During the next quarter, use will be 
made of this program to get successively more complex responses, which will then be used 
In the further study of stability derivative extraction. 

The programmers working on this problem are L. L Mazzola and K. E. Kavanagh. 

L. L. Mazzola 

Aerophyslcs Research Group 

318 C. THREE DIMENSIONAL AERODYNAMIC LEAD PURSUIT STUDY USING WHIRLWIND I AS A REAL 

TIME COMPUTER 

Coding of the problem Is being carried out In three parts: aerodynamic, 
geometric, and display. The first two of those parts have been completed In the sense 
that they give correct results for certain inputs. Currently, work Is being continued 
on the display routine, trying to image correctly on the scope a given sot of points 
representing an outline of the target. After the desired Images are obtained on the 
scope, the Immediate future plans are to program the computation of points outlining 
the target. After the various parts of the problem are chocked out, the parts will be 
Incorporated Into the whole. 
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Programmers working on this proble 
Instrumentation Laboratory. 


are B. Romberg and J. M. Stark of the MIT 


J. M. Stark 

Instrumentation Laboratory 


322 B. MAXIMUM BUBBLE SIZE 


The calculations are completed and will appear in the Sc. D. thesle of Peter 
Griffith, due to be finished June 1956. Later the results will form part of a technical 
report for DIC Project 6627. 

The radius time curves are wanted to aid In the correlation of boiling heat 
flux and maximum heat flux data. Using part of those computer results, a correlation 
has been completed for the maximum heat flux data for pool boiling of a variety of fluids 
and the correlation Is satisfactory. It remains still to complete a correlation tor 
forced convection, aubcooled boiling. Details of these results will be presented In the 
thesis, 

P. Griffith 

Mechanical Engineering Department 

323 N. ANALYSIS OF CLOUD CHAMBER PHOTOGRAPHS 

During this quarter, a tew additional cases were run, completing the first 
stage of the work Intended for this program. Some time in the near future, the program 
will be enlarged to handle some other types of problems. 

The programmer working on this problem is F. Abrams. 

D. 0. Caldwell 

Physics Department 


327 L. PREDICTION ANALYSIS 

This problem deals with the prediction of certain events when information 
giving rise to the prediction has been garbled by Gaussian noise. The solution Is 
primarily geometric and algebraic, with a statistical analysis of the predictions 
appended. The predictions themselves are plotted on the oscilloscope. 

The program Is In working order. It is being used with alternative methods 
of solution to compare results with previously obtained hand-computed cases. 

L. Peterson 
E. Hutcheson 
Lincoln Laboratory 

332 C. GAME-THEORY OPTIMIZATION OF INTERCEPTION SYSTEM 

The problem described In the previous Whirlwind Summary Report was completed 
during the quarter. Entries for the 46-by-46 payoff matrix of the reduced game were 
computed. A study of the matrix Indicated that the reduced-game solution would consist 
of mixed strategies for both players and that both mixed strategies would contain a 
large number of pure strategies. These Indications led to two conclusion. First, an 
exact solution of the reduced game would require computer tlse on the order of year*. 
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second, becsuse the defender's optional strategy wa» known to be a pure strategy, the 
reduced game was an Inadequate representation of the actual game. 

The decision was made, therefore, to attempt an approximate solution of the 
exact game rather than an exact solution of the approximate (reduced) game. * "«thod 
of successive approximations, called the method of fictitious play, was mlopte 
procedure consisted of choosing an alternating succession of strategies for the enemy 
and the defender, each choice being based on the past strategies used by the opponent. 

The method converged rapidly: after 22 iterations an approximate solution was obtained, 
and the true gome value was determined within 1 per cent. The total computing time for 
this run was 29 minutes. 

A final program, lasting 6 minutes, was run to compute some additional results. 
This run yielded data specifying (a) the optimum response function of the cascaded 
director and interceptor and (b) the rms interceptor acceleration of the optimum system. 

The programmer for this problem was K. Kavanagh. 

0. R. Welti 

Dynamic Analysis and Control Laboratory 

333 A. BOUNDARY -LAYER CHARACTERISTICS OF A STEADY LAMINAR FLOW OF A COMBUSTIBLE 
MIXTURE OVER A HOT SURFACE 

During the past few years, the study of ignition and flame stabilization in a 
steady flow of a combustible mixture by means of a hot surface has received increasing 
attention in the combustion field. If a usual bluff-body flame holder is replaced by a 
hot surface, it has been found experimentally that the combustion intensity can be 
increased at least twenty times. However, a basic study of the problem is as yet not 
available. 

A program, both theoretical and experimental, has been set up recently to study 
the basic characteristics of this problem. Four theses are now involved in the experi¬ 
mental program and three term projects are engaged in the theoretical study. Preliminary 
work shows that it is ossentlal to investigate the boundary-layer characteristics of the 
flow of a combustible mixture over a hot surface. This investigation will permit one to 
understand not only the problem of surface ignition in particular but also the effects 
of reaction rate and transfer of momentum, mass, and heat on combustion processes in 
general. 

Partial differential equations of continuity, momentum, and energy are first 
developed for a laminar boundary layer, and then transformed into a series of ordinary 
differential equations. Solutions are to be obtained with the aid of the Whirlwind I 
computer for specific reactions and temperatures of the hot surface. 

Solutions of the first five sets of the differential equations were obtained 
for two values of the temperature of the hot surface and for a second-order reaction of 
the Arrhenius type. 

A study of the profiles of the velocity, temperature, concentration, and 
reaction rate In the boundary layer next to the hot surface indicates that few more 
solutions of the boundary-layer equations are necessary to study the effect of surface 
temperature on incipient ignition. 
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The programmers working on this problem are Professor T. Y. Toong and 
A. Shashaty of the Mechanical Engineering Department. 

T. Y. Toong 

Mechanical Engineering Department 
334 C. PARAMBTRIC STUDY OP MODAL COUPLING AND DAMPING 

(Airing this quarter, the work has been concentrated on the solution of s 
single vibratory differential equation and the simultaneous solution of both two snd 
three coupled vibratory equations. A satisfactory solution has been obtained for each 
of the above three cases based on one set of Initial conditions. To obtain a solution 

to the problem for a different set of Initial conditions, the program must be rerun with 

a modification tape for the constants. Thus a parametric study at the present time 
involves a series of runs with a new modification tape for each run. 

Puture work on the problesi will entail the alteration of the set-up to produce 

a parametric study of the problem In only a single run instead of a series of runs. 

K. Wetmoro 

Aoroelastlc and Structures Research Lab 
338 C. SELECTION OF AIRCRAFT GENERATOR COOLING SYSTEMS 

With increasing flight speeds it has become Impossible to cool aircraft gen¬ 
erators with unprocessed ram air. Systems of turbines, compressors, and heat exchangers 
must be provided to cool the air to the point where it can remove generator losses without 
requiring excessive generator temperatures. Use of such systems imposes a penalty on 
aircraft performance due to the weight of components, and in addition a much larger 
penalty due to drag created by taking air aboard at high velocity and exhausting it at 
low velocity. Equations exist relating this drag to system parameters, and it is desired 
to find the minimum-drag system for an aircraft flying at Mach 3, 70,000 feet. 

A selected range of discreet values of system variables are given to the computer 
together with a program for calculating drag. Another program causes the computer to try 
all combinations of variables and record the drag caused by each. A selection of minimum 
penalty system may then be made by checking weight variations of low-drag systems. 

Five hundred and sixty calculations of drag are required, each calculation 
requiring approximately 200 numerical computations. 

The original problem as described above has been extended to include optimi¬ 
zation of cooling systems for all Mach numbers up to 4 and all altitudes up to 100,000 
feet, as well as investigation of off-design performance of cooling systems. Two types 
of generators have been considered Instead of the original one. Some generator design 
work is also planned for programming. 

The project is substantially completed with the exception of Macb 4 systems. 

R. Moroney 

Energy Conversion Group 

Electrical Engineering Department 
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339 A. BEAM VIBRATION 

An earlier study of the equation 

+ if*. o 

3 x 4 ^ t 2 

was reported on in "Numerical Treatment of a Fourth Order Parabolic Partial Differential 
Equation" by Stephen H. Crandall, J. Assoc. Comp. Mach., 1,111-118 (1954). An improved 
finite difference recurrence formula has now been developed and it is desired to test 
this in actual computations. 

This problem was used as an exercise for the students of MIT Course 2.215, 
Methods of Engineering Analysis, to illustrate the fundamentals of digital computer 
programming. A program was prepared, run and partially debugged by interested students. 


1 expect to complete the program debugging and make a few runs during the next 


quarter. 


S. H. Crandall 

Mechanical Engineering Department 
341 C. STATISTICAL AND DYNAMIC FORECASTING METHODS 

This project is the start of an extended research program utilizing various 
mathematical procedures for the establishment of empirical equations for weather pre¬ 
diction. 

All of the mathematical procedures to be used cannot be stated at present, us 
many are dependent upon initial results. The initial procedure is the computing of the 
elements of covariance matrices. There are several of these matrices to be computed, 
one for each weather element of Interest. Each of these matrices involves input data of 
approximately 10,000 numbers. 

The method of operating on these covariance matrices is in unpublished form, 
and because of its extent it is not feasible to include it here. However, in the not 
too distant future it will be included in a project Scientific Report. 

IXirlng the current quarter the major effort has been devoted toward the 
systematization of a method for multiple linear prediction. In this method a large number 
of predictors is originally chosen. The matrix of the covariances of the predictors is 
obtained. This matrix is then diagonalized. Those linear combinations of the original 
predictors whose variances are the large elements on the diagonalized matrix are then 
used as new predictors. In this way the total number of predictors is reduced. Whirlwind 
I has been used to obtain the necessary cross-products for the covariance matrix. The raw 
data amounts to 10,500 numbers and the elements for a matrix of the 75th order are the 
results. This programming has been under the direction of Mies Elizabeth A. Kelley. 

One phase of the Project is concerned with an extensive study of the extended- 
range forecasting possibilities of the methods developed by the Synoptic Climatology 
Project, now defunct. Tho sea-level pressure and 700 mb height anomaly charts for each 
March from 1946 through 1955 have been expressed in torras of 26 coefficients of 
Tschebyscheff orthogonal polynomials. All maps cover a 135 point qrid bounded by 25°N 
and 65°N an*J by 45°W and 185°W. The points have been taken at each 10° of longitude and 
each 5° of latitude. The Whirlwind computer has been utilized in order to determine the 
coefficients. This work has been under the supervision of Mi. William D. Sellers. 
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Work is being continued in numerical experiments to investigate the climatic 
statistics of very simplified dynamic models of the atmosphere. The dynamic models are 
defined by a series of non-linear difference equations, which express the chsnge with 
respect to time of the amplitude of certain harmonic components of a flow pattern. 
Integration la carried forward by short time atepe. The programming on this phase of 
the project has been done by Mr. Kirk Bryan. 

All of the work done on the project has been supervised by Professor Kdward 
N. Lorenz, MIT Department of Meteorology. 

E. A. Kelley 
Meteorology Department 

343 C. WEATHER PREDICTION 

The success of the research conducted by the Synoptic Climatology Group at 
MIT (Whirlwind Problem No. 155 N.) has indicated that a statistical approach to the 
weather forecast problem should be fully explored. It Is planned to Incorporate meteor¬ 
ological knowledge and theory Into the statistical procedure. The research results of 
our project have indicated fruitful lines of approach. 

There are two aspects of the problem (1) It Is planned to derive linear 
prediction equations utilizing the circulation parameters developed by Wadsworth and 
Bryan and used by the Synoptic Climatology Group. Theory and observation have demon 
strated that the vertical structure of sea-level pressure systems have a direct Influence 
on their behavior and, therefore, on tho future weather. For this reason a non-linear 
approach to the prediction problem will be made through the classification of back data 
Into various classes of atmospheric structure. Linear prediction equations will be devci 
oped for each class of situation Instead of attempting to derive a universal 1 Inoar pre¬ 
diction equation as has been attempted In the past. 

The 1Inoar prediction equations are 
14 14 


% + £ V.! t 1 

1=1 1 1,1 1=1 


Vo.l + c ?I + "o * eS ! + ,8 0 


where £ is the weather element, Z's, P's and S's 
are the coefficients to be obtained. 


are the known variables and a' 


The values of E, Z's S's and P's are already computed. Whirlwind will solve 
for the parameters a Q , u's-f. 

The method of obtaining the coefficients is Least 8quares Analysis of the 
Multiple Linear Regression Equation. The solution of the resulting matrix Is by the 
method described by Crout. It Is possible to obtain both the coefficients and significant 
statistical parameters from this method. 

Once tho coefficients are obtained and the resulting statistical parameter ur*. 
studied, the equation will be evaluated substituting variables independent of these used 
in determining the coeffleients. These results will be verified with observed weather 

conditions. 
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(2) The second aspect of the problem Is the development of new predictors to 
replace the Z's In the above equation. These predictors will be so chosen to represent 
the maximum amount of variance of a field with a minimum number of parameters. Whirlwind 
will perform standard arithmetic operations to obtain the latent roots and latent vectors 
of symmetric matricos. The elements of these matrices will be the correlation coefficients 
or covariances between the Z's. 

During the past quarter, some preliminary work on a desk calculator was under¬ 
taken to classify the data for aspect (1) of the problem. Plans are now under way to 
make an initial computation on this non-linear approach to prediction. 

A study was undertaken to determine sets of orthogonal functions to describe 
patterns of atmospheric variables. Essentially, the technique requires diagonallzation 
of a covariance matrix of such variables. Use was made of Whirlwind library subroutines 
in the initial programming for this problem. 


Since a minimum is desired for the sum of all discounted costs and since an optimal 
(minimum cost) production rate will be chosen at each period, we have. 
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The last equation la found by virtue of the firat-order difference relation: 


The future work on prediction will utilize these new orthogonal functions. 
T .e programmers working on this problem are E. Kelley and B. Shorr. 
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J. M. Austin 
Meteorology Department 


which says that the difference in inventory between the (n+l)st and the nth period is 
the excess of nth period production over nth period demand. 


344 B. AN INVENTORY AND PRODUCTION CONTROL MODEL 

Initial study has been centered in solving numerically an inventory and 
production model where three distinct quadratic (by quadratic we mean the form ax^+bx 4 »c) 
costs are considered: (1) The cost of production levels, (2) The cost of inventory 
levols, and (3) The cost of changes ln the production rate. 

If we have, 
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* the minimum cost of operating for a given number of 
periods when one starts the nth period with inventory, 
x n , and the production rate of the previous period is 

•n-1* 

* demand of nth period; a random variable 

« expected value of ( ) with respect to y n 

* cost during the nth period of production and inventories. 
In general, this function will be a sum of the costs of 
type (1), (2), and (3) mentioned above. 

s discount rate 


and, if we were to study the business for only one period, a minimum cost production 
rate for the nth period could be found by minimizing the expected value of g with respect 
to * . Hence. n 
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Tb® right hand sldo of the equation ia, of 
In business Is that one of considering the 


course. Independent of A common practise 

effect of future periods, n+l, n+2, and so on. 


In the event that f Q conalders the expected costs of many future periods and 
that g n Is quadratic ln form, one finds an analytic expression tor f n such that. 
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A x 2 + B r 2 + C x * , 
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+ D x + S x , + F, 
n n-1 


where the solutions to A, B, C, D, E, F, are given by six non-linear equations of the 
fora, 


A 2 + f <B,...F) A + g 1 (B,...F) = 0 
B 2 + f 2 <A,C, .. .F) B + g 2 (A,B, .. .F) = 0 


F 2 + f g (A,...E) F + g g (A,B,...E) - 0 
The f's and g's are non-linear functions of the indicated variables. 

Whirlwind I programs are being used for the calculation of the coat. f n , and 
the optimal decision for production, z n , for the quadratic model we have outlined. One 
hopes that the answers to this first phase of the problem will allow us to mak# mors 
general assumptions about the cost function, g n . An obvious case is the nscesslty of non 
symmetrical costs for positive and negative inventories as well as costs which are moro 
disproportionately expensive than quadratic costs. An example is the region of high 
production rate changes where one feels a need for an exponential-type coat. 


R. M. Oliver 

Operations Research Group 


32 


S3 


¥ 

i 


I 







T 


APPROVED FOR PUBLIC RELEASE. CASE 06-1104. 


WHIRLWIND CODING AND APPLICATIONS 


345 B. MATRIX MULTIPLICATION 

The problem of the buckling of a clamped shallow spherical thin elastic shell 
subject to a uniform external pressure can bo studied by means of the boundary value 
problem. (See references 1. and 2.) 

4 _ < t .. + a r . 4 »+«. px + 

k 2 x X 2 

<g" + JB' - 4 e) -t - -1/2 t 2 /x 
k x 

x - 0: t, g regular 


1 : 


0, g' - ng 


where k Is d geometric parameter, and P a load parameter. A convenient method for the 
solution o 1 the problem involves assuming t and g in the form of series of complete 
orthogonal functions and by moans of a conversion of the non-linear problem to a sequence 
of linear problems to get the coefficients In the series. This leads to sets of equations 
of the form, 


-<r ) 


t (p) + b (p) = r + c iJ a (p ; 1) b (p ; 1) 
n n n n 1 J 


,‘ p > + (Ji ,2 b ‘P> 
n k n 


M (P-1) (P-1) 


The matrix calculations in the right hand member can be conveniently programmed and 
used over again and again at each stage of the perturbation and also for different values 
of the parameter k that are of interest. 

Solutions have been obtained over a certain range of the parameter k during 
this quarter. It is planned to complete the range of k corresponding to important 

physical applications. 

A. Archor 

Mechanical Engineering Department 

Hoterences 

1. Rolssnor, E., Proc. of Symp. In Appl. Math., Amor. Math. Soc., Vol. 3 (1950) 

2. Simons, R, U., Doctoral Thesis MIT (Math), 1955. 

347 B. SOLUTION OP SETS OP TEN LINEAR SIMULTANEOUS EQUATIONS 

The problem Involves the solution of twelve sets of linear algebraic 
simultaneous equations by the Crout Reduction method. 

Bach of these sots of equations will solve for the ultimate settlements of 
tho footings of a rigid structure, (which In this case are the settlements of the 
footings of the John Dorranco Laboratory) for a soil having a particular compressibility 
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and a superstructure having a particular rigidity. 


The object of my theeis Is to investigate a structure of different rigidities 
and a soil of different compressibilities and solve for the differential settlements for 
each of these combinations and the secondary stresses which these settlements will cause 
In the members of the structure. 

This problem was successfully completed In the first run. 

S. G. Fattal 

Civil Engineering Department 

348 A. UNSTEADY SUPERSONIC PLOWS WITH SPHERICAL SYMMETRY 

The equations which govern the spherically symmetric motion of a compressible 

fluid are: 


A + ( /* u) r + t l° u * 0 
<^>u) t + (P +^u 2 ) r + 2 ^>U 2 o 0 

(+ \ + pu +/ ^r \ + f ( ft pu + ) - 0 

where /o , u, and p are the density, velocity, and pressure respectively at a point 
distant r from the centre of the gas at time t; subscripts denote differentiation. 

The equations 1, 2 and 3 are replaced by finite difference equations by the 
following substitutions: 


( 1 ) 

( 2 ) 

(3) 


3f t 

*t ‘ 't 


» r 


*P ~ 2 (f R + V 

S3 t 


2 ( g R ft S > 


At 


jA r 


Ar 


where f , f , f are the values of f at point P, R and S. 

P R 3 

This form of the forward derivative (f*) Introduces an artificial viscosity 
term in the difference equations so that discontinuities In the variables (which may 
arise due to the presence of shocks) are smoothed over. 

The resulting system of three finite difference equations, was solved numer¬ 
ically using Whirlwind I. g 

Starting from an Initial distribution (at t«0) of /« «1 + Sw * r ,p«(l+5e ) 
u = 0, values of /* , p and u were found on a rectangular mesh of A r » .05 .At-.0125 
from r=0 to 4 and t*0 to 3. 

This requires a total of 19,200 repetitions of the basic cycle. 

During the last quarter, the problem ha. been coded and preliminary result, have 
been obtained for two types of initial conditions, 
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Numerical results show that the Initial 
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-4r 2 -4r 2 

a = 1 + 5e and » 1 + 2e , at t * 

outgoing wav© leavos a region of overexpansion at the centre which causes an ingoing 
shock wave to occur. This shock wav© is reflected from the centre as an outgoing wave 
which finally subsides. 

During the coming quarter it is intended to develop a method whereby certain 
errors, introduced by the numerical procedure, may be eliminated. 

L. Roberts 

Mathematics Department 

349. SOLUTION OP PARTIAL DIPPERENTIAL EQUATIONS 

This problem is concerned with the solution of a set of 23 non-linear partial 
differential equations which arose in the study of a transient heat transfer problem. 
Pirst-order finite difference approximation equations were used to solve this problem. 

Since the basic program contained 1773 instructions, it was handled by the 
core memory, and a solution time of 29 seconds per time level was realized. Since there 
were 400 time levels, a single solution required about 196 minutes (3.3 hours). 

As part of this study, this same problem was solved on the IBM Type 6S0 Magnetic 
Drum Calculator. The main 650 program contained 1743 instructions and solution for a 
single time level was obtained in 29 minutes. The 650 program was not optimum programmed, 
if it had been, the solution time would have been reduced to seven minutes for a single 
time level. 

The results obtained on Whirlwind and the Type 650 Calculator agree to six 
decimal digits--which is quite good considering the basic differences In the machine 
structure, etc. 

Programmers working on this problem were Dr. P. M. Verzuh and M. Hermann of 
the Office of Statistical Services, and A. Zabludowsky of the Digital Computer Laboratory. 

P. M. Verzuh 

Office of Statistical Services 

351 B. NON-UN1PORM PUEL DISTRIBUTION 
Non-Uniform Puel Loading in R eactors 

Reactors loaded in a uniform pattern throughout the core allow only a partial 
exploit of the outer core regions, as maximum heat generation at core center Is limited 
by heat transfer and metallurgical properties. 

A deliberate redistribution of fuel, moderator, neutron absorbers and coolant 
will modify the neutron flux and equalize maximum temperatures and stresses within the 

reactor. 

Tiie Age-Diffusion equations are solved for 1 energy groups, N core-zones and 
M reflector-sones. Material constants vary for each group i and each region N. The 
Age-Diffusion equations are approximated by finite differences. The T^placlan In the 
leakage term becomes a second derivative when introducing modified fluxes pi . r 0 1 
in spherical geometry. 
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Boundary condition* impose P» to vanish at the center ot the core and at the 
extrapolated boundary. Neutron currents must be equal at core reflector Interface By 
successive Iterations the distribution of fuel concentration In the core produclng’an 
arbitrary thermal flux can bo approximated. 

Age-Diffusion equations in finite differences: (solved by Vhlrlwlnd I) 


F* . - k 1 F 1 + F 1 + I 1 
trt-1 n n n-1 n 


(n = 0,.,.NfM) 


The first throe terms includo leakage and absorption. Tile last term describes the slowing 
down and production of neutrons. Tho number of Iterations is a function of the accepted 
error. 

The problem Is divided Into two major programs: space criterion and criticality 

criterion. 

Space Criterion 

From an arbitrary assumed thermal flux the correct spatial flux distribution 
is approximated by successive Iterations (regardless of criticality conditions). 

The program for the space criterion Is being fulfilled for a 2-group, 2-reglon 
reactor with 18 core zones and 20 reflector zones. Print-out of single iterations shows 
rapid convergence, A program for any number of boundaries will be developed. 

Criticality Criterion 

Material properties of the reactor are changed to acquire prescribed criticality. 
The changes are fuel concentration, fuel distribution, reflector thickness, reflector 
composition, and critical size. 

The routine for criticality criterion will be programmed according to the 
various changes In material properties as mentioned above. 

The programmers working on this problem are A. Sutter and U. Troost. 

A. Sutter 

Nuclear Engineering 

352 B. CRITICAL WHIRLING FREQUENCIES OF PROPELLER SHAFTS 

Tho whirling or lateral vibrations of a propeller shafting Induced by the 
blade frequency of the propeller will be calculated by numerical Iteration. The method 
Involves a numerical solution to a fourth degree differential equation. Frequencies 
will be assumed and the residual moment at the free end of the shafting will be plotted 
versus these frequencies. The point at which the end moment equals zero Is a resonant 
point or natural whirling frequency of the shafting Bystem. This method of output was 
chosen to enable comparison with existing data. The numerical calculations are: 


- 1 g s 0 


y > C x sinAx + C 2 cos\x + C 3 slnh)jc + C 4 coshkx 
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Using end conditions to solve for constants the following parameters are 
obtained. These constants may be considered lumped parameters. 


w„ l; B, I, 

p P011 

3 3 

*1 L x E o x o 


TT 4 ( - 2 - ) 4 
1 ' 180 ' 




J L 0 E 1 *1 ( I 1 8 } 

* '180^ 


"l L 1 S 0 X o 


“l X 1 


dl 


W L 3 
1 1 


<ii|- >’ f VTc 


al 


Vl 


01 


dl 


cosfl slnhfl - sing coshp 
slnhg - sing 

cosg CQBhg - 1 
gCsinhg - sing) 


2g sing sinhg 
slnhg - sing 


B„ I 
0 0 




fu, V 


V 


N+l 


0 .. 


N.W-1 


, . Cf . 

3 J 18 

[<Vn Vi,n + c .n <Vn Vi.n] 
fo o , "n-1.N ( Vn ) 

t-1 Vl.N + -j 


ftirlng the past quarter a program has been written for this problem and is in 
the process of being tested. 

Programmers working on this problem ore Lt. C. K. Drandt, USN, Lt. J. C. Snyder, 
USN, and Lt. C. K. Thompson, USC0, associated with the HIT Department of Naval Architecture. 

C. K. Drandt 
Naval Architecture 
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353 A. WAITING LINK8 - CONSTANT HOLDING TINS 
The problem la aa follow** 

n -x 

Given: x « »u, a(n,x) « 5-2-— 

n; 


n « 0, 1. 


to compute: 11m P(n,i) for n » 0, 1, 

1 -* °° 


where 


P(n,0) • 1 


for all n and P(n,l+1) » £ P(srt-J,i) a(n-J) 

J-0 


It is known that for u< 1 the limits exist. Computations are required for m - 1, 2, S, 
10, 20 and u = .2, .4, .6, .7, .8, .85, .90, .95, .98, and .99. 

The limiting probabilities here are the steady state probabilities for the 
uumber of units waiting or being served In a queue with m servicing channels, Poisson 
arrivals, and constant holding (servicing) time, u is the ratio of mean arrival rate to 
total mean servicing rate. 

The feasibility has been verified and seven problems have been successfully 


In the future, lt is planned to run up to 43 remaining problems. 

H. P. Galllher 
Operations Research Group 
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APPENDIX 


APPENDIX 


1. SYSTEMS ENGINEEKING 


Performance Records for the Whirlwind Computer System 

The WWI computer system has been performing a 168 hour week since 1 June 1955. 
This 7 day per week, 24 hours per day schedule Is expected to be continued for at least 
one year. The amelioration of computer system reliability noted In the two most recent 
quarterly periods was apparently a direct result of the refinements made to several of 
the equipment areas during the preceding 6 months. 

The following table illustrates how the equipment has performed during the 
most recent 18 months. 

w 

28 Sept.'54 7 Oct.'55 30 Dec.'55 

to to to 

6 Oct.'55_29 Dec.'55_8 Mar.'56 


Total Computer Operating Time 
In Hours 

Total Time Lost In Hours 

Percentage Operating Time Usable 

Average Uninterrupted Operating 

Time Between Pailure Incidents 
In Hours 

Total Number of Failure Incidents 

Failure Incidents per 24-Hour Day 

Average Lost Time Per Incident 
In Minutes 

Average Preventive Maintenance Time 
Per Day In HourB 


7537 1876 1584 


230 

48 

36.3 

97.0 

97.4 

97.7 

12 

17.8 

22.2 

605 

103 

70 

1.92 

1.32 

1.06 

36.6 

28 

22.6 

1.55 

1.7 

1.45 









3. VISITORS 

Tours of the *hlrlwlnd I Installation Include a shoving of the film, "Making 
Electrons Count", a computer demonstration, and an Informal discussion of the major 
computer components, ftirlng the past quarter, the following 10 groups totalling 143 
people visited the computer Installation: 

January 4 
January 6 
January 10 
January 11 
January 16 
January 19 
January 30 
January 30 
March 7 

March 15 

The procedure of holding Open House at the Digital Computer Laboratory on the 
first and third Tuesday of each month has continued, with the exception of March 20th 
when the Institute was closed because of a severe blizzard. A total of 112 people 
attended the five Open House tours during the quarter, representing members and friends 
of the MIT Community, Children's Medical Center, Medfield State Hospital, and Harvard 
University. 


MIT Aeroelastlc Laboratory 

MIT Course In Fundamentals In Electrical Engineering 

Lexington High School 

MIT Course In "Pulsed-Data Systems 

Harvard Business School 

MIT Air Force Cadets 

General Electric Co. 

"Frontiers of Science" Group 

Dr. Heinz Krekeler, German Ambassador to the United States, 
and members of his staff 
Mills College (New York) 


During the last quarter, there were also 45 Individuals who made brief tours 
of the computer installation at different times. Represented by the individuals were: 
MIT, Melpar, Lynn Mutual Insurance Co., A.O.Smith Corp., Barnard College, International 
Business Machines Corp., the United states Navy, Burroughs Corporation, Princeton 
University, Metal Box Co., Ltd., Continental Can Co., Harvard University, Chance Vought 
Aircraft, Tohoka University, Oberlln College, and Babson Institute. 
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MACHINE METHODS OF COMPUTATION AND NUMERICAL ANALYSIS 
Faculty Supervisors : 


Philip M. Morae, Chairman 
Samuel H. Caldwell 
Herman Feshbach 
Jay w. Forreater 
James B. Keswick 
Chla Chlao Lin 


R esearc h Associate: 


Bayard Rankin 

Research Assistan ts: 

Joseph Hershenov 
M. Douglas Mcllroy 
James W. Schleslnger 
Hubertus J. Welnltschke 
Fernando J. Corbato 
Harvey Fields 
Zoltan Fried 
Raymond F. Stora 
Arnold Tubis 
Aaron Temkin 
Jack L. Uretsky 
Jukka A. Lehtinen 
Joseph B. Walsh 
Marius Troost 
Norman F. Ness 
Ralph G. Gray 


PROJECT whirlwind 


Physics 

Electrical Engineering 
Physics 

Electrical Engineering 
Mechanical Engineering 
Mathematics 


Mathematics • 


Mathematics 

Mathematics 

Mathematics 

Mathematics 

Physics 

Physics 

Physics 

Physics 

Physics 

Physics 

Physics 

Mechanical Engineering 
Mechanical Engineering 
Chemical Engineering 
Geology 

Civil Engineering 


Digital Computer^Laboratory" Cl - ntlflC peering Computation. Group at the 

Frank M. Verzuh, Head 
Dean N. Arden 
Rheldon P. Best (Abs.) 

John M. Frankovlch 
Frank C. Helwlg 
leonard Roberts 
Jack Koseman 
Arnold Siegel 
Murray Watkins 
Monroe R. Weinstein 
Alexander Zabludowsky 
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